THIS CONTAINS A PROGRAM AND DATA FOR USING THE MMC2 MODEL IN APPLICATION TO
CLUSTERS OF HYDROGEN MOLECULES.
THIS FILE CONTAINS 4 ITEMS:
1) A DATA FILE WHICH IS READ AS FORTRAN LOGICAL FILE NUMBER 195
2) A DATA FILE WHICH IS READ AS FORTRAN LOGICAL FILE NUMBER 196
3) A DATA FILE WHICH IS READ AS FORTRAN LOGICAL FILE NUMBER 197
4) A FORTRAN PROGRAM FOR MODEL POTENTIAL EVALUATION.
THE PROGRAM IS EXPLAINED BY THE COMMENT LINES WITHIN. IT MAY BE NECESSARY TO
READ THROUGH THE WHOLE THING TO UNDERSTAND HOW TO USE IT. IT CAN BE COMPILED,
LINKED, AND EXECUTED ON ITS OWN AND IT WILL RUN A TEST CALCULATION. THE VALUE
IT SHOULD PRODUCE IS GIVEN IN COMMENT LINES. IF THE CORRECT VALUES ARE NOT
PRODUCED, IT MEANS THERE IS SOME FORTRAN COMPILER OPTION THAT IS REQUIRED.
(THE PROGRAM IS ESSENTIALLY WRITTEN AS FORTRAN77 CODE, BUT THERE ARE
VARIATIONS IN FORTRAN77 COMPILERS AND OPTIONS FROM ONE SYSTEM TO ANOTHER.)
THE MAIN PROGRAM CAN BE REMOVED AND THE REMAINING SUBROUTINES USED BY
SOMEONE'S OWN PROGRAM TO EVALUATE H2-CLUSTER ENERGIES. THE INFORMATION THAT
NEEDS TO BE PASSED TO THE SUBROUTINES IS INFORMATION ABOUT THE CLUSTER
GEOMETRY. THAT MEANS THE (X,Y,Z) COORDINATES OF THE H2 CENTERS IN ATOMIC
UNITS, AND THEIR EULER ANGLES OF ORIENTATION WITH RESPECT TO A FIXED
LABORATORY AXIS SYSTEM (WHICH CAN BE DEFINED BY ONE OF THE MOLECULES, IF
DESIRED).
BREAK THIS FILE INTO SEPARATE DATA FILES AND THEN THE PROGRAM FILE. EACH FILE
SHOULD START WITH THE LINE FOLLOWING THE BLANK LINE THAT IDENTIFIES THE
ITEM.
ITEM 1 - DATA FILE 195 ........... FILE STARTS AFTER BLANK LINE
3 THIS VALUE (FORMAT-I5) IS THE NUMBER OF CENTERS HEREIN.
HH CENTER OF MASS
1.651 419.03 0.0 0.0
HYDROGEN IN HH
0.0 151.525 1.00783 0.0
EXTRA CENTERS IN HH
0.879 0.0 0.0 0.0
ITEM 2 - DATA FILE 196 ........... FILE STARTS AFTER BLANK LINE
1 :: NUMBER OF MOLECULES HEREIN-FORMAT(I5). (X,Y,Z) IN A.U.'S
HH
5
1 0.0 0.0 0.0 1
2 -0.702033204 0.0 0.0 0
2 0.702033204 0.0 0.0 0
3 -1.079978402 0.0 0.0 0
3 1.079978402 0.0 0.0 0
ITEM 3 - DATA FILE 197 ........... FILE STARTS AFTER BLANK LINE
HH CENTER OF MASS
0.0 0.0 0.0 0.0
-0.536221043 0.0 0.0
0.0 -0.782250254 0.0
0.0 0.0 -0.7822502540
6.472765753 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 4.644939136 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 4.644939136
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0 0 0 3.867560549 0 0 0 0.662362243 0 0 0 0.662362243
0 0 0 0 1.420154381 0 1.420154381 0 0 0 0 0
0 0 0 0 0 1.420154381 0 0 0 1.420154381 0 0
0 0 0 0 1.420154381 0 1.420154381 0 0 0 0 0
0 0 0 0.662362243 0 0 0 2.582574613 0 0 0 1.049758077
0 0 0 0 0 0 0 0 0.766400805 0 0.766400805 0
0 0 0 0 0 1.420154381 0 0 0 1.420154381 0 0
0 0 0 0 0 0 0 0 0.766400805 0 0.766400805 0
0 0 0 0.662362243 0 0 0 1.049759077 0 0 0 2.582574613
HYDROGEN IN HH
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
EXTRA CENTERS IN HH
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
ITEM 4 - FORTRAN PROGRAM ............ FILE STARTS AFTER BLANK LINE
PROGRAM TEMPDRIVER
IMPLICIT REAL*8 (A-H,O-Z)
COMMON / CNTRLM3 / NTYPES(100), NCENTERS(100), IATOM(30,30),
1 IATYPE(30,30), STRUCTURE(30,30,3)
DIMENSION U(27), GEOM(3,6)
C
C THIS PROGRAM ROUTINE IS DESIGNED TO DEMONSTRATE HOW TO USE
C THE SET OF SUBROUTINES THAT FOLLOW. THIS ROUTINE "TEMPDRIVER"
C IS DESIGNED TO RUN AN EXAMPLE. IT IS INTENDED THAT IT BE
C REMOVED AND THE REMAINING SUBROUTINES LINKED WITH ANY OTHER
C PROGRAM, PROVIDED COMMON BLOCKS ARE ARRANGED, AND CALLING
C PARAMETERS TO THE ROUTINE "MMCPOT" ARE SET UP.
C
C SOME CARE MUST BE TAKEN FOR THE COMPATIBILITY OF ONE'S
C FORTRAN COMPILER WITH THE USE OF COMMON BLOCKS HEREIN.
C DATA IS COMMUNICATED FROM SUBROUTINE TO SUBROUTINE VIA
C THE COMMON BLOCKS.
C
C THE SUBROUTINE "MAKEMOLS" NEEDS TO BE CALLED ONCE TO INPUT
C PARAMETERS AND PROPERTIES. THE FOLLOWING STATEMENT SHOULD BE
C INCLUDED IN ANY PROGRAM THAT TAKES THE PLACE OF "TEMPDRIVER".
C
CALL MAKEMOLS
C
C >>> EXAMPLE OF DATA TO FEED TO THE "MMCPOT" ROUTINE VIA COMMON
C BLOCKS ETC. (SET UP IN THE ROUTINE THAT TAKES THE PLACE OF
C "TEMPDRIVER")
C
C 1. THREE MOLECULES (3 H2'S FOR THIS EXAMPLE)
C
NMOLECS = 3
C
C 2. SPECIFY WHICH TYPES OF MOLECULES (IDENTICAL IN THIS
C DEMONSTRATION)
C
NTYPES(1) = 1
NTYPES(2) = 1
NTYPES(3) = 1
C
C 3. SET UP THE "U" MATRICES (A 3 X 3 FOR EACH MOLECULE). FOR THIS
C EXAMPLE THEY WILL BE SET UP TO BE THE IDENTITY (I.E., NO
C REORIENTATION; EULER ANGLES = 0.0) FOR THE FIRST AND THE THIRD.
C THE SECOND MOLECULE WILL HAS A U-MATRIX CORRESPONDING TO A
C 45-DEGREE ROTATION ABOUT THE Z-AXIS
C
DO 10 I=1,27
U(I) = 0.0D 0
10 CONTINUE
U(1) = 1.0D 0
U(5) = 1.0D 0
U(9) = 1.0D 0
RR = 1.0D 0 / SQRT(2.0D 0)
U(10) = RR
U(11) = -RR
U(13) = RR
U(14) = RR
U(18) = 1.0D 0
U(19) = 1.0D 0
U(23) = 1.0D 0
U(27) = 1.0D 0
C
C 4. POSITION OF THE THREE MOLECULES WILL BE AT (0, 0, 0),
C (6.6, 0, 0), AND (0.0, 6.6, 0) [ATOMIC UNITS]
C SPECIFYING EULER ANGLES IS NOT NECESSARY FOR THIS
C EXAMPLE (ONLY), SINCE THAT INFORMATION
C IS CONVEYED VIA THE "U" MATRICES
C
GEOM(1,1) = 0.0D 0
GEOM(1,2) = 0.0D 0
GEOM(1,3) = 0.0D 0
GEOM(2,1) = 6.6D 0
GEOM(2,2) = 0.0D 0
GEOM(2,3) = 0.0D 0
GEOM(3,1) = 0.0D 0
GEOM(3,2) = 6.6D 0
GEOM(3,3) = 0.0D 0
C
C>>> END OF "DATA"
C
C>>> OUPUT: THIS RUN SHOULD PRODUCE THE FOLLOWING
C
C THE MMC2 INTERACTION ENERGY IS -0.87429438E-04 A.U.
C -19.188576 CM-1
C>>>
C
CALL MMCPOT(VVVV,NMOLECS,U,GEOM)
VVCM = VVVV * 219475.0D 0
WRITE(6,900) VVVV, VVCM
900 FORMAT(' THE MMC2 INTERACTION ENERGY IS',G26.8,' A.U.'/
1 ' ',G26.8,' CM-1')
END
SUBROUTINE MMCPOT(VVVV,NMOLECS,U,GEOM)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION U(3,3,NMOLECS), GEOM(NMOLECS,6)
C
C THIS ROUTINE IS AN ENCODING OF AN MMC-STYLE INTERACTION POTENTIAL,
C WHICH IS DESIGNATED MMC2.
C A CALL TO THIS SUBROUTINE YIELDS A VALUE FOR "VVVV" WHICH IS THE
C POTENTIAL ENERGY (IN A.U.) FOR THE GIVEN ARRANGEMENT OF MOLECULES.
C
C THE ARRAY "U" IS A SET OF TRANSFORMATION MATRICES, ONE FOR EACH OF
C THE NMOLECS (NUMBER OF MOLECULES). THESE ARE TRANSFORMATION
C MATRICES THAT ROTATE THE MOLECULE FROM ITS REFERENCE ORIENTATION
C TO ITS CURRENT ORIENTATION IN SPACE.
C
C "GEOM" HOLDS TRADITIONAL X, Y, Z, AND EULER ANGLE (R, S, T)
C COORDINATES FOR EACH MOLECULE. THIS ROUTINE USES ONLY THE
C X, Y, Z COORDINATES TAKING THE NECESSARY ORIENTATION INFORMATION
C FROM THE U-MATRICES.
C
C VERSION 1.0 OCTOBER 10, 1995
C VERSION 1.1 OCTOBER 19, 1995
C VERSION 1.2 JULY 17, 1997 (DISTRIBUTION II)
C VERSION 1.3 AUGUST 20, 2003 (ADAPTED FOR H2-CLUSTERS)
C
COMMON / CNTRLF / EFIELD(3)
COMMON / CNTRLM1 / PMOMS(13,200), POLAR(144,200)
COMMON / CNTRLM2 / C6(200), C12(200), AMASS(200)
C
C THE FIRST DIMENSION FOR THE ARRAYS IN /CNTRLM1/ IS THE MAXIMUM
C NUMBER OF TYPES OF ELECTRICAL CENTERS. IN /CNTRLM2/ IT IS THE
C MAXIMUM NUMBER OF 6-12 CENTERS.
C
C THE ORDERING OF THE FIRST INDEX IN "PMOMS" (PERMANENT MOMENTS) IS
C (1) CHARGE
C (2), (3), (4) DIPOLE MOMENT
C (5) ... (13) QUADRUPOLE MOMENT
C THE ORDERING IN "POLAR" (POLARIZABILITIES) IS SUCH THAT THE
C 144 ELEMENTS FOR EACH CENTER TYPE AS IF CAST AS A
C 12-BY-12 ARRAY, THE ROWS AND COLUMNS WOULD BE THOSE OF
C ELEMENTS 2-13 OF "PMOMS".
C
COMMON / CNTRLM3 / NTYPES(100), NCENTERS(30), IATOM(30,30),
1 IATYPE(30,30), STRUCTURE(30,30,3)
C
C "NTYPES" IS DIMENSIONED TO BE THE MAXIMUM NUMBER OF INTERACTING
C SPECIES. THE VALUE OF NTYPE(I) IS THE TYPE-NUMBER FOR THE I-TH
C MOLECULE.
C "NCENTERS" IS DIMENSIONED TO BE THE MAXIMUM NUMBER OF MOLECULE
C TYPES (30) CONSTRUCTED AND USED IN A GIVEN RUN. THE VALUE OF
C NCENTERS(I) IS THE NUMBER OF CENTERS IN THE I-TH TYPE OF MOLECULE.
C "IATOM" IS DIMENSIONED AS THE NUMBER OF MOLECULE TYPES BY THE
C MAXIMUM NUMBER OF ATOMS IN A MOLECULE. IATOM(I,J) IS THE VALUE
C ASSOCIATED WITH THE J-TH ATOM (CENTER) IN THE I-TH MOLECULE
C TYPE, AND ITS VALUE IS
C 0 IF IT IS A 6-12 CENTER (ONLY)
C 1 IF IT IS AN ELECTRICAL/6-12 CENTER
C 2 IF IT IS AN ELECTRICAL CENTER (ONLY)
C "IATYPE" IS LIKE "IATOM" BUT THE VALUES ARE THE CENTER TYPES
C USED WITH PMONS, POLAR, C6 AND C12.
C "STRUCTURE" HOLDS INFORMATION ON THE STRUCTURE OF EACH MOLECULE
C TYPE. IT IS DIMENSIONED THE SAME AS "IATOM" BUT TAKES ON 3
C VALUES CORRESPONDING TO THE X,Y,Z POSITION COORDINATES OF THE
C J-TH ATOM IN THE I-TH MOLECULE TYPE, RELATIVE TO THE MASS CENTER
C OF THAT MOLECULE TYPE.
C
DIMENSION POSITION(100,30,3), XXX(3), VATK(13)
DIMENSION QL(100,30), DL(3,100,30), TL(3,3,100,30)
DIMENSION UT(3,3), VTEMP(13)
C
C "POSITION" HOLDS THE ROTATED X,Y,Z POSITION COORDINATES OF THE
C CENTERS (SECOND INDEX) IN EACH OF THE MOLECULES (FIRST INDEX).
C "XXX" IS THE VECTOR FROM THE L-CENTER TO THE K-CENTER IN THE LOOPS.
C "VATK" IS THE POTENTIAL POLYTENSOR AT THE K-CENTER.
C
DO 100 K=1,NMOLECS
KTYPE = NTYPES(K)
KCENT = NCENTERS(KTYPE)
DO 80 KK=1,KCENT
KATYPE = IATYPE(KTYPE,KK)
CALL TRANSF1(U(1,1,K),PMOMS(2,KATYPE),DL(1,K,KK))
CALL TRANSF2(U(1,1,K),PMOMS(5,KATYPE),TL(1,1,K,KK))
DO 60 I=1,3
SS = 0.0D 0
DO 40 J=1,3
SS = SS + U(I,J,K) * STRUCTURE(KTYPE,KK,J)
40 CONTINUE
POSITION(K,KK,I) = SS + GEOM(K,I)
60 CONTINUE
80 CONTINUE
100 CONTINUE
C
C NEXT IS A "DOUBLE LOOP" TO GO OVER ALL PAIRS OF CENTERS ON
C DIFFERENT MOLECULES. THE ELECTRICAL POTENTIAL AND 6-12
C POTENTIALS ARE EVALUATED TOGETHER.
C
V66 = 0.0D 0
V12 = 0.0D 0
VEL = 0.0D 0
DO 800 K=1,NMOLECS
KTYPE = NTYPES(K)
KCENT = NCENTERS(KTYPE)
CALL TRANSPOSE(U(1,1,K),UT)
DO 780 KK=1,KCENT
KATOM = IATOM(KTYPE,KK)
KATYPE = IATYPE(KTYPE,KK)
DO 210 M=1,13
VATK(M) = 0.0D 0
210 CONTINUE
DO 700 L=1,NMOLECS
IF(K .EQ. L) GO TO 700
LTYPE = NTYPES(L)
LCENT = NCENTERS(LTYPE)
DO 680 LL=1,LCENT
LATOM = IATOM(LTYPE,LL)
LATYPE = IATYPE(LTYPE,LL)
R2 = 0.0D 0
DO 240 M=1,3
XXX(M) = POSITION(K,KK,M) - POSITION(L,LL,M)
R2 = R2 + XXX(M) * XXX(M)
240 CONTINUE
R2 = 1.0D 0 / R2
IF( KATOM .EQ. 2 .OR. LATOM .EQ. 2) GO TO 300
IF(L .GE. K) GO TO 300
C
C EVALUATE 6-12 TERM
C
R6 = R2 ** 3
R12 = R6 * R6
V66 = V66 - C6(KATYPE) * C6(LATYPE) * R6
V12 = V12 + C12(KATYPE) * C12 (LATYPE) * R12
300 CONTINUE
IF(KATOM .EQ. 0 .OR. LATOM .EQ. 0) GO TO 680
R1 = SQRT(R2)
R3 = R1 * R2
R5 = R3 * R2
C
C USE ROTATED PERMANENT MOMENTS (ONLY) OF THE CENTERS IN MOLECULE-L
C FROM "QL": THE CHARGE ON THE LL-CENTER.
C FROM "DL": THE DIPOLE MOMENT ON THE LL-CENTER.
C FROM "TL": THE QUADRUPOLE MOMENT ON THE LL-CENTER
C
CALL FIELD(XXX,R1,R2,R3,R5,VATK,QL(L,LL),
1 DL(1,L,LL),TL(1,1,L,LL))
680 CONTINUE
700 CONTINUE
IF(KATOM .EQ. 0) GO TO 780
C
C ADD THE EXTERNAL FIELD THAT IS SPECIFIED IN /CNTRLF/
C
VATK(2) = VATK(2) + EFIELD(1)
VATK(3) = VATK(3) + EFIELD(2)
VATK(4) = VATK(4) + EFIELD(3)
C
C USE THE TRANSPOSE OF THE K-MOLECULE'S TRANSFORMATION MATRIX
C TO TRANSFORM THE FIELD AND GRADIENT >>BACK<< TO
C THE MOLECULE FRAME SO AS TO AVOID TRANSFORMING THE POLARIZABILITY
C TENSORS.
C CALL "POLENER" TO EVALUATE THE ENERGY OF POLARIZATION.
C
CALL TRANSF1(UT,VATK(2),VTEMP(2))
CALL TRANSF2(UT,VATK(5),VTEMP(5))
CALL POLENER(VELK,VTEMP,POLAR(1,KATYPE),PMOMS(1,KATYPE))
VEL = VEL + VELK
780 CONTINUE
800 CONTINUE
VVVV = VEL + V66 + V12
RETURN
END
SUBROUTINE FIELD(XXX,R1,R2,R3,R5,VV,Q,D,T)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XXX(3), D(3), T(3,3), VV(13), C(3), DD(3)
C
C FOR A VECTOR "XXX" DRAWN FROM SOME ELECTRICAL POINT (WITH
C CHARGE "Q", DIPOLE "D" AND QUADRUPOLE "T"), THIS ROUTINE
C FINDS THE FIELD AND FIELD GRADIENT AT THE END OF THE VECTOR.
C THE FIELD IS STORED AS ELEMENTS 2,3,4 OF "VV" AND THE
C FIELD GRADIENT IS STORED AS ELEMENTS 5 TO 13.
C
C R1, R2, R3, R5 ARE INVERSE POWERS OF THE LENGTH OF THE VECTOR.
C
RD = 0.0D 0
RTR = 0.0D 0
TRT = 0.0D 0
DO 30 I=1,3
S = 0.0D 0
RD = RD + XXX(I) * D(I)
TRT = TRT + T(I,I)
DO 20 J=1,3
S = S + T(I,J) * XXX(J)
20 CONTINUE
RTR = RTR + XXX(I) * S
C(I) = S * R2
DD(I) = D(I) + 10.0D 0 * C(I)
30 CONTINUE
A = R2 * (RD - TRT)
B = R2 * R2 * RTR
A3B = A + 3.0D 0 * B
A5B = 3.0D 0 * (A + 5.0D 0 * B)
A7B = 5.0D 0 * (A + 7.0D 0 * B)
QA5B = Q + A5B
QA7B = Q + A7B
R53 = 3.0D 0 * R5
R3Q = R3 * (Q + A5B)
VV(1) = VV(1) + R1 * (Q + A3B)
NN = 1
DO 60 I=1,3
VV(I+1) = VV(I+1) + R3 * ( D(I) + 6.0D 0 * C(I)
1 - XXX(I) * QA5B )
NN = NN + 3
VV(NN+I) = VV(NN+I) - R3Q
DO 50 J=I,3
N = NN + J
VV(N) = VV(N) + R53 *
1 ( XXX(I) * XXX(J) * QA7B
2 - XXX(I) * DD(J)
3 - XXX(J) * DD(I) + 2.0D 0 * T(I,J) )
50 CONTINUE
60 CONTINUE
VV(8) = VV(6)
VV(11) = VV(7)
VV(12) = VV(10)
RETURN
END
SUBROUTINE MAKEMOLS
IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE READS IN THREE DATA FIELDS THAT ULTIMATELY SERVE
C TO CONSTRUCT THE MOLECULE TYPES OUT OF INDIVIDUAL CENTERS.
C THE INDIVIDUAL CENTERS MAY BE ON ATOMS BUT MAY BE ANYWHERE
C ELSE AS WELL.
C
COMMON / CNTRLM1 / PMOMS(13,200), POLAR(144,200)
COMMON / CNTRLM2 / C6(200), C12(200), AMASS(200)
COMMON / CNTRLM3 / NTYPES(100), NCENTERS(30), IATOM(30,30),
1 IATYPE(30,30), STRUCTURE(30,30,3)
DIMENSION TITLE(80), TITLE2(80)
C
C "NCURRENT" IS THE CURRENT NUMBER OF CENTERS FOR WHICH THERE ARE
C 6-12 AND/OR ELECTRICAL PARAMETERS.
C
REWIND 195
REWIND 196
REWIND 197
READ(195,800) NCURRENT
DO 100 N=1,NCURRENT
READ(195,900) TITLE
READ(197,900) TITLE2
DO 10 I=1,80
IF(TITLE(I) .EQ. TITLE2(I)) GO TO 10
WRITE(6,910) N, TITLE, TITLE2
GO TO 20
10 CONTINUE
20 READ(195,*) C6(N), C12(N), AMASS(N), AXT
C
C "AXT" IS NOT USED IN THIS VERSION.
C
READ(197,*) (PMOMS(I,N),I=1,13)
READ(197,*) (POLAR(I,N),I=1,144)
100 CONTINUE
C
C "NUMTYPS" IS THE NUMBER OF MOLECULE TYPES REPRESENTED ON
C FILE 196. THE MAXIMUM IS 10. THE POTENTIAL WILL BE EVALUATED
C FOR A LARGE COLLECTION OF MOLECULES, EACH OF WHICH IS ONE OF
C THE "NUMTYPS" OF MOLECULES DESCRIBED VIA FILE 196.
C
READ(196,800) NUMTYPS
DO 300 N=1,NUMTYPS
READ(196,900) TITLE
READ(196,*) NUMCENT
NCENTERS(N) = NUMCENT
DO 200 K=1,NUMCENT
READ(196,*) IATYPE(N,K), (STRUCTURE(N,K,I),I=1,3),
1 IATOM(N,K)
200 CONTINUE
300 CONTINUE
C
800 FORMAT(I5)
900 FORMAT(80A1)
910 FORMAT(' >>> WARNING: FOR DATA SET',I4,' LABELS ON FOR195'
1 /' AND FOR197 DO NOT MATCH'/5X,80A1/5X,80A1)
END
SUBROUTINE POLENER(ENER,VFIELD,POLAR,PMOMS)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION VFIELD(13), PMOMS(13), POLAR(12,12)
C
C THIS ROUTINE EVALUATES THE ADDITIVE ELECTROSTATIC ENERGY OF SOME
C CENTER GIVEN THE PERMANENT MOMENTS ("PMOMS"), THE POLARIZABILITIES
C ("POLAR") AND THE POLYTENSOR FIELD ("VFIELD").
C
ENER = VFIELD(1) * PMOMS(1) * 2.0D 0
DO 20 I=2,13
II = I - 1
X = PMOMS(I)
DO 10 J=2,13
JJ = J - 1
X = X - POLAR(II,JJ) * VFIELD(J)
10 CONTINUE
ENER = ENER + X * VFIELD(I)
20 CONTINUE
ENER = 0.5D 0 * ENER
RETURN
END
SUBROUTINE TRANSF1(U,V,Y)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION U(3,3), V(3), Y(3)
DO 20 I=1,3
S = 0.0D 0
DO 10 J=1,3
S = S + U(I,J) * V(J)
10 CONTINUE
Y(I) = S
20 CONTINUE
RETURN
END
SUBROUTINE TRANSF2(U,A,B)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION U(3,3), A(3,3), B(3,3), X(3,3)
DO 30 I=1,3
DO 20 J=1,3
S = 0.0D 0
DO 10 K=1,3
S = S + U(I,K) * A(K,J)
10 CONTINUE
X(I,J) = S
20 CONTINUE
30 CONTINUE
DO 60 I=1,3
DO 50 J=1,3
S = 0.0D 0
DO 40 K=1,3
S = S + X(I,K) * U(J,K)
40 CONTINUE
B(I,J) = S
50 CONTINUE
60 CONTINUE
RETURN
END
SUBROUTINE TRANSPOSE(U,UT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION U(9), UT(9)
UT(1) = U(1)
UT(2) = U(4)
UT(3) = U(7)
UT(4) = U(2)
UT(5) = U(5)
UT(6) = U(8)
UT(7) = U(3)
UT(8) = U(6)
UT(9) = U(9)
RETURN
END