Program ROTENER C********************************************************************** C* Release September 1996 * C* * C* Program ROTENER * C* * C* by M. Nardelli * C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica * C* e Chimica Fisica della Universita' di Parma, Centro di Studio per * C* la Strutturistica Diffrattometrica del CNR. * C* Viale delle Scienze 78, I-43100 Parma, Italy * C* E-mail: nardelli@ipruniv.cce.unipr.it * C* * C* Literature: J.Chem.Soc.Dalton, 1990, p.179; ibid. 1991, p.203 * C* Gazz. Chim. Ital., (1991), 121, 433-439 * C* QCPE BULLETIN. (1991), 11, 61 * C* Bull. Chem. Soc. Jpn.(1994). 67, 189-195 * C* J.Chem.Inform. Computer Sci. (1994). 34, 972-975 * C* Acta Cryst.(1994). C50, 1323-1326 * C* * C* This program calculates the difference molecular potential energy * C* after rotations of one or two molecular fragments about given * C* directions. The first two atoms of each set to be rotated define * C* the directions about which the rotations occur. * C* Two options are available: isolated molecule or molecule in the * C* crystal. In the latter case the crystal environment is considered * C* static. * C* The van der Waals energy is calculated only if the atoms involved * C* in the contacts are: H, C, N, O, F, Si, P, S, Cl, Fe, Br, Sn, I, * C* Pb. * C* The Coulombic energy is considered only if the atomic partial * C* charges are given with non zero values. * C* The options isolated molecule and molecule in the crystal (static * C* environment are available. * C* * C* INPUT * C* 1. TITLE (A80) * C* 2. Bravais symbol (A1) * C* 3. KY code for atomic species. KY=0 if not given, KY=1 if given * C* 4. a, b, c, alpha, beta, gamma (free FORMAT) * C* 5. NNS = number of atomic species (free FORMAT) * C* 6. Atomic species between apices (free FORMAT) * C* 7. Format of atomic coordinates (fractional) and charges * C* (electrons) [e.g.(1X,6A1,1X,3F9.5,F7.3)] * C* 8. N (total number of atoms) * C* 9. Atom, x, y, z, q (6A1,3F8.5,F7.3) (max number of atoms: 400) * C*10. The number of the atoms of the first rotating fragment,theta * C* minimum (rotation angle, degr), theta maximum, delta theta * C* (step of theta)(free format). * C*11. The labels of four atoms defining a torsion angle about the * C* first direction of rotation (4(6A1)) * C*12. The labels of the atoms of the first rotating fragment. The * C* origin of the reference orthogonal system of axes is on the * C* first atom and the rotation occurs counterclock-wise about the * C* direction defined by the first two atoms (12(6A1)/12(6A1)) * C* N.B. If the number of atoms is EQUAL to 12 or 24..., put a blank * C* record after their list * * C*13. The number of the atoms of the second rotating fragment,theta * C* minimum (rotation angle, degr), theta maximum, delta theta * C* (step of theta) (free format). * C*14. The labels of four atoms defining a torsion angle about the * C* second direction of rotation (4(6A1)). * C*15. The labels of the atoms of the second rotating fragment. The * C* rotation occurs counterclock-wise about the direction defined * C* by the first two atoms (12(6A1)/12(6A1)). * C*16. Data for environment in the crystal (not given if isolated * C* molecule is considered: IC,NE,NT (Free Format) * C* IC=1 if the space group is acentric * C* =-1 if the space group is centric * C* NE=Number of the equivalent positions (the general, x y z, the * C* centrosymmetric and those generated by the Bravais lattice * C* must be omitted * C* NT=Maximum translation required (suggested 1 when the set of * C* atoms in the general position is the nearest to the origin, * C* 2 in the other cases) * C*17. Equivalent position: one equivalent position per record (20A1) * C* Examples: 1/2-X,1/2+Y,1/2+Z (not 0.5-X, etc...) * C* Y-X,-X,1/3+Z (not 0.3333+Z) * C* * C* N.B.- If the second fragment is part of the firts one the names of * C* the atoms of the latter fragment must be present also among * C* those for the first one. * C* * C* The program interactively ask for the number of fragments to * C* be considered and for isolated molecule or in the crystal. * C* For rotations of 360 degr. and one fragment, steps of 5 degr. * C* are suggested, while with two fragments, steps of 10 degr. are * C* advisable. Of course the computing time depends also upon the * C* number of atoms. * C* * C* OUTPUT * C* * C* In user defined files: * C* * C* Title * C* * C* On request: crystal data and atomic coordinates. * C* * C* * C* In the case of one rotating fragment: * C* ===================================== * C* * C* A table giving the rotation angles, the torsion angles, the * C* difference energy with respect the unrotated molecule, in kcal/mol * C* and kJ/mol. In addition, interatomic contacts less than a given * C* value (requested interactively; suggested 2-2.5 A) and energies * C* corresponding to each rotation of the fragment are also quoted. * C* * C* On request: tables of data for two-dimensional graphic routines, * C* the names of the files are given by the user. * C* * C* * C* In the case of two fragments: * C* ============================ * C* * C* A table giving the difference energy (with respect the unrotated * C* molecule) for each couple of rotation angles. * C* * C* A line-printer plot of the difference energy values as function of * C* the rotation angles. * C* * C* On request: a table giving the T1, T2, DE data suitable as input * C* for a tridimensional plotting routine. * C* * C* ----------------- * C* WARNING * C* GETTIM and GETDAT are intrinsec routines of the IBM Professional * C* FORTRAN Compiler PROFORT Version 1.0 [(C) Copyright IBM Corp * C* 1984, (C) Copyright Ryan-McFarland Corp 1984] library used to * C* get time and date. If a different compiler is used, the * C* following sentences of the Main must be changed or deleted. * C* At the beginning: * C* C-----Read the initial time * C* CALL GETTIM(IHR,IMN,ISC,IHS) * C* HS=IHS * C* TM0=(IHR*3600)+(IMN*60)+ISC+HS/100.0 * C* At the end: * C* C-----Print computing time and date * C* CALL GETDAT(IYR,IMT,IDY) * C* CALL GETTIM(IHR,IMN,ISC,IHS) * C* HS=IHS * C* TM=(IHR*3600)+(IMN*60)+ISC+HS/100.0 * C* IF(TM.LT.TM0)TM=TM+86400 * C* DTM=TM-TM0 * C* WRITE(IO,'(/'' Total computing time: '',F9.2,'' s'')')DTM * C* WRITE(IO,140) * C* WRITE(IO,140)IDY,MONTH(IMT),IYR,IHR,IMN,ISC,IHS * C* 140 FORMAT(I4,'/',A4,'/',I4,3X,I4,':',I2,':',I2,'.',I2/) * C* * C********************************************************************** CHARACTER*1 TIT(80),FRT(80),AT(500,6),ATR1(500,6),MINUS(2),LR, 1 AT1(4,6),AT2(4,6),ATR2(500,6),LP,QS,QL(2),BRV,AAA1(500,6), 2 AAA2(500,6),AFR1(500,6),AFR2(500,6),STPE(38,20),SDL1(2,500), 3 LA1(2,20),AAAA(6),LBA(2),CTEMP(6) CHARACTER*14 FILGRF,FILVIS,FILSRF CHARACTER*4 MONTH(12) CHARACTER*2 SDL2(500),LA(20) INTEGER*2 IHR,IMN,ISC,IHS,IYR,IMT,IDY DIMENSION XRO1(500,3),XRO2(500,3),XB(500,3),QB(500), 1 R1(4,3),TAU1(180),TAU2(180),ENRG(72,72),IO1(72),IO2(72), 2 IENRG(72,72),XYZ(5184,3),R2(4,3),RR1(3),RR2(3),XO1(500,3), 3 XO2(500,3),XO3(500,3),XOO1(500,3),XOO2(500,3),X1K(3), 4 TH(500),DEJ1(500),NT1(4),NT2(4),XF1(500,3),XXO(500,3), 5 XBO(500,3),NA4(500),QF(500) DOUBLE PRECISION BI1(3,3),BI2(3,3) EQUIVALENCE (LA(1),LA1(1,1)),(SDL2(1),SDL1(1,1)) COMMON/GRP1/O(3,3) 1 /GRP2/AT 2 /GRP3/MC(400,8) 3 /GRP5/XO(500,3) 4 /GRP6/NA1(500),NA2(500),NA3(500) 5 /GRP7/Q1(500),Q2(500),Q3(500) 6 /GRP8/DMAX 7 /GRP10/X(500,3),Q(500) 8 /LBL1/C(3),A(3) 9 /LBL3/IC,NE,NT COMMON/LBL4/STPE 2 /INPO/IN,IO 3 /GRG1/SDL2 4 /GRG2/LA DATA QL/'N','n'/ DATA MONTH/'Jan.','Feb.','Mar.','Apr.','May.','Jun.','Jul.', 1 'Aug.','Sep.','Oct.','Nov.','Dec.'/ C-----Read the initial time CALL GETTIM(IHR,IMN,ISC,IHS) HS=IHS TM0=(IHR*3600)+(IMN*60)+ISC+HS/100.0 C-----Constants and I/O files RD=0.0174532925 IN=13 IO=15 CALL INPOUT C-----Read title, codes and parameters READ(IN,100)TIT 100 FORMAT(80A1) WRITE(IO,101)TIT 101 FORMAT(1X,80A1) READ(IN,'(A1)')BRV READ(IN,*)KY READ(IN,*)C,A CALL ORTHOG(C,A) DO 60 I=1,20 60 LA(I)=' ' IF(KY.EQ.0)GO TO 59 READ(IN,*)NNS READ(IN,*)(LA(I),I=1,NNS) DO 57 I=1,NNS MINUS(1)=LA1(1,I) MINUS(2)=LA1(2,I) CALL MINUMA(MINUS) LA1(1,I)=MINUS(1) 57 LA1(2,I)=MINUS(2) 59 READ(IN,100) FRT WRITE(*,102) 102 FORMAT(' Do you want crystal data and atomic coordinates', 1 ' printed? (Y/N)'/' ?>') READ(*,'(A1)')LP IF(LP.EQ.QL(1).OR.LP.EQ.QL(2))GO TO 1 WRITE(IO,103)(C(I),A(I),I=1,3) 103 FORMAT(' Crystal data'/' a = ',F7.4,25X,'alpha= ',F6.2/ 1 ' b = ',F7.4,25X,'beta = ',F6.2/' c = ',F7.4,25X,'gamma= ',F6.2) WRITE(IO,104)((O(I,J),J=1,3),I=1,3) 104 FORMAT(/' Orthogonalization matrix'/3(1X,3F10.5/)) WRITE(IO,105) 105 FORMAT(' Atomic coordinates'//' Atom',6X,'X/a',6X,'Y/b', 1 6X,'Z/c',12X,'X',8X,'Y',8X,'Z'/) 1 READ(IN,*)N DO 4 I=1,N READ(IN,FRT)(AT(I,L),L=1,6),(X(I,J),J=1,3),Q(I) DO 2 K=1,5,2 MINUS(1)=AT(I,K) MINUS(2)=AT(I,K+1) CALL MINUMA(MINUS) AT(I,K)=MINUS(1) 2 AT(I,K+1)=MINUS(2) C-----Orthogonal coordinates DO 3 J=1,3 XO(I,J)=0.0 DO 3 K=1,3 3 XO(I,J)=XO(I,J)+O(J,K)*X(I,K) C-----Matrix of the atom symbols DO 58 L=1,6 58 AAAA(L)=AT(I,L) CALL LABL(AAAA,KY,NNS,LA1,LBA,KKD) SDL1(1,I)=LBA(1) SDL1(2,I)=LBA(2) IF(LP.EQ.QL(1).OR.LP.EQ.QL(2))GO TO 4 WRITE(IO,106)(AT(I,L),L=1,6),(X(I,J),J=1,3),(XO(I,K),K=1,3) 106 FORMAT(1X,6A1,2X,3F9.5,5X,3F9.4) 4 CONTINUE C-----Set up the connectivity matrix MC(400,8) CALL BONDI(N) C-----Read in the rotation parameters WRITE(*,107) 107 FORMAT(' Key in the number of fragments to be rotated (1 or 2)'/ 1 ' ?>') READ(*,*)KC WRITE(*,108) 108 FORMAT(' Key in 0 if the molecule is considered isolated '/ 1 ' 1 if the molecule is considered in the crystal'/ 2 ' ?>') READ(*,*)KDE CALL READAT(NG1,T1,T3,DT1,AT1,ATR1) IF(KC.EQ.2) CALL READAT(NG2,T2,T4,DT2,AT2,ATR2) C-----Consider crystal environment IF(KDE.EQ.1)THEN READ(IN,*)IC,NE,NT DO 5 I=1,NE 5 READ(IN,'(20A1)')(STPE(I,L),L=1,20) NEE=NE C-----Environment of the first rotating fragment NR1=NG1-2 DO 6 I=3,NG1 DO 6 L=1,6 6 AFR1(I-2,L)=ATR1(I,L) CALL CONTACT(BRV,N,NR1,AFR1,NL,AAA1,XXO,QF) C-----Environment of the second rotating fragment IF(KC.EQ.2)THEN NR2=NG2-2 DO 7 I=3,NG2 DO 7 L=1,6 7 AFR2(I-2,L)=ATR2(I,L) NE=NEE CALL CONTACT(BRV,N,NR2,AFR2,NB,AAA2,XBO,QB) ENDIF ENDIF T1=T1*RD T3=T3*RD DT1=DT1*RD IF(KC.EQ.2)THEN T2=T2*RD T20=T2 T4=T4*RD DT2=DT2*RD ENDIF C-----First rotating fragment DO 11 I=1,NG1 DO 10 J=1,N DO 8 L=1,6 IF(ATR1(I,L).NE.AT(J,L)) GO TO 10 8 CONTINUE DO 9 K=1,3 XO1(I,K)=XO(J,K) 9 XF1(I,K)=X(J,K) NA1(I)=J Q1(I)=Q(J) GO TO 11 10 CONTINUE WRITE(*,109)(ATR1(I,L),L=1,6) 109 FORMAT(' Atom ',6A1,' does not exist in the molecule.'/ 1 ' Check your data. Calculation must stop.') STOP 11 CONTINUE IF(KC.EQ.2)THEN C-----Second rotating fragment DO 15 I=1,NG2 DO 14 J=1,N DO 12 L=1,6 IF(ATR2(I,L).NE.AT(J,L)) GO TO 14 12 CONTINUE DO 13 K=1,3 13 XO2(I,K)=XO(J,K) NA2(I)=J Q2(I)=Q(J) GO TO 15 14 CONTINUE WRITE(*,109)(ATR2(I,L),L=1,6) STOP 15 CONTINUE ENDIF C-----Unrotating fragment KL=0 DO 19 J=1,N DO 16 I=3,NG1 IF(NA1(I).EQ.J)GO TO 19 16 CONTINUE IF(KC.EQ.2)THEN DO 17 I=3,NG2 IF(NA2(I).EQ.J)GO TO 19 17 CONTINUE ENDIF KL=KL+1 NA3(KL)=J Q3(KL)=Q(J) DO 18 K=1,3 18 XO3(KL,K)=XO(J,K) 19 CONTINUE NG3=KL C-----Atoms involved in the torsion angles CALL ATORS(N,AT1,NG1,NT1,NA1,R1) IF(KC.EQ.2)CALL ATORS(N,AT2,NG2,NT2,NA2,R2) C-----Add crystal environment to the fix fragment IF(KDE.EQ.1)THEN CALL ADD(N,NL,AAA1,XO3,XXO,Q3,QF,NG3) C-----Add environmental atoms to the atom symbol matrix DO 62 I=1,N DO 61 L=1,6 61 CTEMP(L)=AT(I,L) CALL LABL(CTEMP,KY,NNS,LA1,LBA,KKD) SDL1(1,I)=LBA(1) 62 SDL1(2,I)=LBA(2) IF(KC.EQ.2)THEN NI=0 DO 23 I=1,NB DO 23 J=1,NL DO 20 L=1,3 IF(ABS(XBO(I,L)-XXO(J,L)).GT.0.0001)GO TO 23 20 CONTINUE NI=NI+1 DO 22 K=I,NB-NI DO 21 L=1,3 XBO(K,L)=XBO((K+1),L) QB(K)=QB(K+1) 21 AAA2(K,L)=AAA2((K+1),L) DO 22 L=4,6 22 AAA2(K,L)=AAA2((K+1),L) 23 CONTINUE NB=NB-NI CALL ADD(N,NB,AAA2,XO3,XBO,Q3,QB,NG3) C-----Add environmental atoms to the atom symbol matrix IF(N.GT.500) WRITE(*,'(/'' WARNING''/'' The total number '', 1 ''of atoms is greater than 500''/)') DO 64 I=1,N DO 63 L=1,6 63 CTEMP(L)=AT(I,L) CALL LABL(CTEMP,KY,NNS,LA1,LBA,KKD) SDL1(1,I)=LBA(1) 64 SDL1(2,I)=LBA(2) ENDIF ENDIF C-----Calculation of the molecular energy for the unrotated molecule KD=0 CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E00) KD=1 CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E01) IF(KC.EQ.2)THEN KD=2 CALL ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,E02) E0=E00+E01+E02 ELSE E0=E00+E01 ENDIF E0J=E0*4.184 IF(KDE.EQ.0)WRITE(IO,110)E0,E0J 110 FORMAT(/' Non-bonded energy for the isolated molecule:'/ 1 ' E0 = ',F10.2,' kcal/mol'/' E0 = ',F10.2,' kJ/mol') C-----Ask for a upper limit value of delta(E) WRITE(*,111) 111 FORMAT(' Key in the maximum value of delta(E), kJ/mol'/' ?>') READ(*,*)DEMAX IF(KC.EQ.1)THEN WRITE(*,112) 112 FORMAT(/' Key in the maximum value of the distance (A)', 1 ' to have the contacts quoted'/' ?>') READ(*,*)DMAX C-----Write the head of the output WRITE(IO,113) 113 FORMAT(/' Rotation of the fragment') WRITE(IO,114)((ATR1(I,L),L=1,6),I=1,NG1) 114 FORMAT(12(1X,6A1)/12(1X,6A1)) WRITE(IO,115)((ATR1(I,L),L=1,6),I=1,2),((AT1(I,L),L=1,6),I=1,4) 115 FORMAT(' about the ',6A1,'-',6A1,' direction [Phi(1)]'/ 1 ' Torsion ',3(6A1,'-'),6A1,2X,'[Tau(1)]'/) WRITE(IO,116) 116 FORMAT(' Rotations are counter-clockwise'/) WRITE(IO,117)((AT1(I,L),L=1,6),I=1,4) 117 FORMAT(' Non-bonded difference molecular energy'//' Rotation',7X, 1 'Torsion',16X,'DE'/' angle',8X,2(6A1,'-'), 2 5X,'kcal/mol',5X,'kJ/mol'/15X,2('-',6A1)/) ELSE WRITE(IO,113) WRITE(IO,114)((ATR1(I,L),L=1,6),I=1,NG1) WRITE(IO,115)((ATR1(I,L),L=1,6),I=1,2),((AT1(I,L),L=1,6),I=1,4) WRITE(IO,113) WRITE(IO,114)((ATR2(I,L),L=1,6),I=1,NG2) WRITE(IO,118)((ATR2(I,L),L=1,6),I=1,2),((AT2(I,L),L=1,6),I=1,4) 118 FORMAT(' about the ',6A1,'-',6A1,' direction [Phi(2)]'/ 1 ' Torsion ',3(6A1,'-'),6A1,2X,'[Tau(2)]'/) WRITE(IO,116) ENDIF C-----Ask for the complete list of data IF(KC.EQ.2)THEN WRITE(*,119) 119 FORMAT(' Do you want the complete list of data? (Y/N)'/' ?>') READ(*,'(A1)')LP IF(LP.EQ.QL(1).OR.LP.EQ.QL(2)) GO TO 24 C-----Write the head of the table giving the complete list of the data WRITE(IO,120) 120 FORMAT(/' Rotation angles',5X,'Torsion angles',9X, 1 'Difference energy'/2X,'Phi(1)',2X,'Phi(2)',5X,'Tau(1)',2X, 2 'Tau(2)',10X,'kcal/mol kJ/mol') C-----Ask for the data to use with graphics 24 WRITE(*,121) 121 FORMAT(' Do you want the list of T1, T2, DE data '/ 1 ' to use as input with plotting routines? (Y/N)'/' ?>') READ(*,'(A1)')QS ENDIF CALL NEWCOO(NG1,XO1,BI1,XOO1,RR1) NE1=0 NN1=0 NX=0 C-----Start the loop of the first fragment rotations 25 NE1=NE1+1 CALL ROTCOO(NG1,XOO1,T1,BI1,XRO1,RR1) IF(KC.EQ.2)THEN C-----Modify the coordinates of the atoms of the second fragment, C-----if this fragment is part of the first one DO 30 I=1,NG1 DO 29 J=3,NG2 DO 26 L=1,6 IF(ATR2(J,L).NE.ATR1(I,L)) GO TO 29 26 CONTINUE DO 28 K=1,3 X1K(K)=0.0 DO 27 KK=1,3 27 X1K(K)=X1K(K)+BI1(K,KK)*XOO1(I,KK) 28 XO2(J,K)=X1K(K)+RR1(K) GO TO 30 29 CONTINUE 30 CONTINUE ENDIF C-----The torsion angle for the first fragment is calculated CALL FRTORS(NG1,NT1,NA1,AT1,NN1,TAU1,XRO1,R1) KD=1 CALL ENERMOL(KC,KD,NG1,NG2,NG3,XRO1,XRO2,XO3,E1) IF(KC.EQ.1) THEN TH(NN1)=T1/RD EJ=E1*4.184 E01J=E01*4.184 DE=E1-E01 DEJ1(NN1)=EJ-E01J IF(DEJ1(NN1).GT.DEMAX)DEJ1(NN1)=DEMAX DE=DEJ1(NN1)/4.184 WRITE(IO,122)TH(NN1),TAU1(NN1),DE,DEJ1(NN1) 122 FORMAT(F7.1,7X,F8.1,3X,2F12.2) WRITE(16,123)TH(NN1) 123 FORMAT(' Rotation angle:',F7.1,' degr.'/) ELSE NE2=0 NN2=0 T2=T20 CALL NEWCOO(NG2,XO2,BI2,XOO2,RR2) ENDIF C-----Start the loop of the second fragment rotations IF(KC.EQ.2)THEN 31 NE2=NE2+1 CALL ROTCOO(NG2,XOO2,T2,BI2,XRO2,RR2) CALL FRTORS(NG2,NT2,NA2,AT2,NN2,TAU2,XRO2,R2) C-----Calculate the non bonded energy for fragment 2 KD=2 CALL ENERMOL(KC,KD,NG1,NG2,NG3,XRO1,XRO2,XO3,E2) DE=(E1+E2)-(E01+E02) DEJ=DE*4.184 IF(DEJ.GT.DEMAX)DEJ=DEMAX DE=DEJ/4.184 TH1=T1/RD TH2=T2/RD IF(LP.EQ.QL(1).OR.LP.EQ.QL(2)) GO TO 32 C-----Printing of the complete table of energy calculations WRITE(IO,124)TH1,TH2,TAU1(NN1),TAU2(NN2),DE,DEJ 124 FORMAT(2(F7.1,1X),3X,2F8.1,6X,2F10.1) 32 ENRG(NE1,NE2)=DEJ NX=NX+1 XYZ(NX,1)=TH1 XYZ(NX,2)=TH2 XYZ(NX,3)=DEJ C-----Rounding off the angular values IF(IO1(NE1).GE.0)IO1(NE1)=TH1+0.5 IF(IO1(NE1).LT.0)IO1(NE1)=TH1-0.5 IF(IO2(NE2).GE.0)IO2(NE2)=TH2+0.5 IF(IO2(NE2).LT.0)IO2(NE2)=TH2-0.5 T2=T2+DT2 IF(T2.GT.T4) GO TO 33 GO TO 31 C-----End of the second fragment rotation loop 33 ENDIF T1=T1+DT1 IF(T1.GT.T3)GO TO 34 GO TO 25 C-----End of the first fragment rotation loop 34 IF(KC.EQ.2)THEN WRITE(IO,'(1X)') C-----Scaling of the energy data EMX=ENRG(1,1) EMN=ENRG(1,1) DO 35 I=1,NE1 DO 35 J=1,NE2 IF(ENRG(I,J).GT.EMX)EMX=ENRG(I,J) IF(ENRG(I,J).LT.EMN)EMN=ENRG(I,J) 35 CONTINUE EINT=EMX-EMN FT=1 IF(ABS(EINT).LE.10)FT=10 IF(ABS(EINT).GT.100)FT=0.1 IF(ABS(EINT).GT.1000)FT=0.01 IF(ABS(EINT).GT.10000)FT=0.001 IF(ABS(EINT).GT.100000)FT=0.0001 WRITE(IO,125)FT 125 FORMAT(' The following energy values have been multiplied by', 1 F8.4) C-----Rounding off the energy values DO 36 I=1,NE1 DO 36 J=1,NE2 IF(ENRG(I,J).GE.0.0)IENRG(I,J)=ENRG(I,J)*FT+0.5 IF(ENRG(I,J).LT.0.0)IENRG(I,J)=ENRG(I,J)*FT-0.5 36 CONTINUE C-----Printing of the three-dimensional energy table IF(NE2.GT.20)GO TO 38 WRITE(IO,126)(IO2(I),I=1,NE2) 126 FORMAT(/' phi2³',20I4) WRITE(IO,127) 127 FORMAT(' _____³',80('_')/' phi1 ³') WRITE(IO,128) 128 FORMAT(6X,'³') DO 37 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,NE2) 129 FORMAT(I5,1X,'³',20I4) 37 WRITE(IO,128) GO TO 50 38 IF(NE2.GT.40)GO TO 41 WRITE(IO,126)(IO2(I),I=1,20) WRITE(IO,127) WRITE(IO,128) DO 39 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20) 39 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=21,NE2) WRITE(IO,127) WRITE(IO,128) DO 40 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,NE2) 40 WRITE(IO,128) GO TO 50 41 IF(NE2.GT.60)GO TO 45 WRITE(IO,126)(IO2(I),I=1,20) WRITE(IO,127) WRITE(IO,128) DO 42 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20) 42 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=21,40) WRITE(IO,127) WRITE(IO,128) DO 43 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,40) 43 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=41,NE2) WRITE(IO,127) WRITE(IO,128) DO 44 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=41,NE2) 44 WRITE(IO,128) 45 WRITE(IO,126)(IO2(I),I=1,20) WRITE(IO,127) WRITE(IO,128) DO 46 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=1,20) 46 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=21,40) WRITE(IO,127) WRITE(IO,128) DO 47 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=21,40) 47 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=41,60) WRITE(IO,127) WRITE(IO,128) DO 48 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=41,60) 48 WRITE(IO,128) WRITE(IO,126)(IO2(I),I=61,NE2) WRITE(IO,127) WRITE(IO,128) DO 49 I=1,NE1 WRITE(IO,129)IO1(I),(IENRG(I,K),K=61,NE2) 49 WRITE(IO,128) 50 IF(QS.EQ.QL(1).OR.QS.EQ.QL(2))GO TO 51 WRITE(*,130) 130 FORMAT(' Key in the name of the file for graphic routines'/' ?>') READ(*,'(A14)')FILSRF OPEN(UNIT=17,FILE=FILSRF,FORM='FORMATTED') REWIND(UNIT=17) WRITE(17,131)((XYZ(I,K),K=1,3),I=1,NX) 131 FORMAT(2F8.1,1X,F15.2) 51 CONTINUE ELSE WRITE(IO,132)DMAX 132 FORMAT(/' Contacts less than',F5.2,' A ','after rotations of the') WRITE(IO,133)((ATR1(I,L),L=1,6),I=1,NG1) 133 FORMAT(12(1X,6A1)/12(1X,6A1)) WRITE(IO,134)((ATR1(I,L),L=1,6),I=1,2) 134 FORMAT(' fragment about the ',6A1,'-',6A1,' direction'/) REWIND(UNIT=16) DO 52 I=1,10000 READ(16,'(80A1)',END=53)TIT 52 WRITE(IO,'(80A1)')TIT 53 WRITE(*,135) 135 FORMAT(' Do you want data suitable for graphic software? (Y/N)'/ 1 ' ?>') READ(*,'(A1)')LR IF(LR.EQ.QL(1).OR.LR.EQ.QL(2))GO TO 56 WRITE(*,136) 136 FORMAT(' Key in:'/' 0 if D(E) is a function of the rotation', 1 ' angle'/' 1 if D(E) is a function of the torsion angle'/ 2 ' ?>') READ(*,*)KLD WRITE(*,137) 137 FORMAT(' Key in the name of the file for GRAFIT'/' ?>') READ(*,'(A14)')FILGRF OPEN(UNIT=17,FILE=FILGRF,FORM='FORMATTED') REWIND(UNIT=17) DO 55 I=1,NN1 IF(KLD.EQ.1)GO TO 54 WRITE(17,138)TH(I),DEJ1(I) 138 FORMAT(F11.1,1X,F15.3) GO TO 55 54 WRITE(17,138)TAU1(I),DEJ1(I) 55 CONTINUE WRITE(*,139) 139 FORMAT(' Key in the name of the file for VISION'/' ?>') READ(*,'(A14)')FILVIS OPEN(UNIT=18,FILE=FILVIS,FORM='FORMATTED') REWIND(UNIT=18) WRITE(18,'(8F10.4)')(DEJ1(I),I=1,NN1) 56 CONTINUE ENDIF C-----Print computing time and date CALL GETDAT(IYR,IMT,IDY) CALL GETTIM(IHR,IMN,ISC,IHS) HS=IHS TM=(IHR*3600)+(IMN*60)+ISC+HS/100.0 IF(TM.LT.TM0)TM=TM+86400 DTM=TM-TM0 WRITE(IO,'(/'' Total computing time: '',F9.2,'' s'')')DTM WRITE(IO,140) WRITE(IO,140)IDY,MONTH(IMT),IYR,IHR,IMN,ISC,IHS 140 FORMAT(I4,'/',A4,'/',I4,3X,I4,':',I2,':',I2,'.',I2/) CLOSE(IN) CLOSE(IO) CLOSE(16) CLOSE(17) CLOSE(18) STOP END SUBROUTINE INPOUT C-----Assign input and output files CHARACTER*14 FILIN COMMON/INPO/IN,IO WRITE(*,100) 100 FORMAT(' Key in the input file'/' ?>') READ(*,'(A14)')FILIN OPEN(UNIT=IN,FILE=FILIN,FORM='FORMATTED') REWIND(UNIT=IN) WRITE(*,101) 101 FORMAT(' Key in the output file'/' ?>') READ(*,'(A14)')FILIN OPEN(UNIT=IO,FILE=FILIN,FORM='FORMATTED') REWIND(UNIT=IO) OPEN(UNIT=16,STATUS='SCRATCH',FORM='FORMATTED') RETURN END SUBROUTINE LABL(AAAA,KY,NS,LA1,LBA,KKD) C-----This routine finds the atom species from its label CHARACTER*1 AAAA(6),LA1(2,20),LBA(2),MINU(26),MAIU(26) COMMON/INPO/IN,IO 1 /MIMA/MINU,MAIU IF(KY.EQ.1)THEN DO 1 J=1,NS IF(AAAA(1).EQ.LA1(1,J).AND.AAAA(2).EQ.LA1(2,J))GO TO 3 1 CONTINUE DO 2 K=1,NS IF(AAAA(1).EQ.LA1(1,K))GO TO 4 2 CONTINUE WRITE(IO,100)AAAA 100 FORMAT('Atom ',6A1,' is not among the given species'/ 1 'The first two characters are assumed') LBA(1)=AAAA(1) LBA(2)=AAAA(2) KKD=1 RETURN 3 LBA(1)=LA1(1,J) LBA(2)=LA1(2,J) RETURN 4 LBA(1)=LA1(1,K) LBA(2)=' ' RETURN ENDIF LBA(1)=AAAA(1) DO 5 I=1,26 IF(AAAA(2).EQ.MAIU(I))GO TO 6 5 CONTINUE LBA(2)=' ' RETURN 6 LBA(2)=AAAA(2) RETURN END SUBROUTINE ATORS(N,AT1,NG1,NT1,NA1,R1) C-----Considers the atoms involved in the torsion angles CHARACTER*1 AT(500,6),AT1(4,6) DIMENSION R1(4,3),XO(500,3),NT1(4),NA1(500) COMMON/GRP2/AT 1 /GRP5/XO DO 5 I=1,4 DO 2 J=1,N DO 1 L=1,6 IF(AT1(I,L).NE.AT(J,L))GO TO 2 1 CONTINUE GO TO 3 2 CONTINUE WRITE(*,100)(AT1(I,L),L=1,6) 100 FORMAT(' Atom ',6A1,' does not exists in the molecule.'/ 1 ' Check your data. Calculation must stop.') STOP 3 DO 4 K=1,3 4 R1(I,K)=XO(J,K) NT1(I)=J 5 CONTINUE DO 6 I=3,NG1 IF(NT1(1).EQ.NA1(I)) GO TO 7 6 CONTINUE RETURN 7 IRIS=NT1(1) NT1(1)=NT1(4) NT1(4)=IRIS IRIS=NT1(2) NT1(2)=NT1(3) NT1(3)=IRIS DO 8 K=1,3 RIS=R1(1,K) R1(1,K)=R1(4,K) R1(4,K)=RIS RIS=R1(2,K) R1(2,K)=R1(3,K) 8 R1(3,K)=RIS RETURN END SUBROUTINE FRTORS(NG1,NT1,NA1,AT1,NN1,TAU1,XRO1,R1) C-----Considers the atoms defining the torsion angle involved in the C-----rotating fragment CHARACTER*1 AT1(4,6) DIMENSION NT1(4),NA1(500),TAU1(180),R1(4,3),XRO1(500,3) DO 1 J=1,NG1 IF(NT1(4).EQ.NA1(J))GO TO 2 1 CONTINUE WRITE(*,100) (AT1(NT1(4),L),L=1,6) 100 FORMAT(' Atom ',6A1,' does not exists in the molecule.'/ 1 ' Check your data. Calculation must stop.') STOP 2 DO 3 K=1,3 3 R1(4,K)=XRO1(J,K) CALL TORS(R1,TAU) NN1=NN1+1 TAU1(NN1)=TAU RETURN END SUBROUTINE ENERMOL(KC,KD,NG1,NG2,NG3,XO1,XO2,XO3,ET) C-----This routine calculates the interatomic contacts and the C-----van der Waals energy CHARACTER*1 ATOM(2,2),AT(500,6) CHARACTER*2 KS(2),SDL2(500) DIMENSION XO1(500,3),XO2(500,3),XO3(500,3),X1(3),X2(3) COMMON/GRP2/AT 1 /GRP3/MC(400,8) 2 /GRP4/ATOM 3 /GRP6/NA1(500),NA2(500),NA3(500) 4 /GRP7/Q1(500),Q2(500),Q3(500) 5 /GRP8/DMAX 6 /GRG1/SDL2 ET=0 IF(KD.EQ.0)THEN DO 5 I=1,NG3-1 DO 5 J=I+1,NG3 DO 2 K=1,8 IF(NA3(J).EQ.MC(NA3(I),K))GO TO 5 DO 1 L=1,8 IF(NA3(J).EQ.MC(MC(NA3(I),K),L))GO TO 5 1 CONTINUE 2 CONTINUE DO 3 K=1,3 X1(K)=XO3(I,K) 3 X2(K)=XO3(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 5 IF(D.GT.10.0)GO TO 5 KS(1)=SDL2(NA3(I)) KS(2)=SDL2(NA3(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q3(I)*Q3(J)/D)*331.979 ET=ET+E+EQ 5 CONTINUE ELSEIF(KD.EQ.1)THEN DO 12 I=1,NG3 DO 12 J=3,NG1 DO 7 K=1,8 IF(NA1(J).EQ.MC(NA3(I),K))GO TO 12 DO 6 L=1,8 IF(NA1(J).EQ.MC(MC(NA3(I),K),L))GO TO 12 6 CONTINUE 7 CONTINUE IF(KC.EQ.2)GO TO 9 DO 8 L=3,NG2 IF(NA1(J).EQ.NA2(L))GO TO 12 8 CONTINUE 9 DO 10 K=1,3 X1(K)=XO3(I,K) 10 X2(K)=XO1(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 12 IF(D.GT.10.0)GO TO 12 KS(1)=SDL2(NA3(I)) KS(2)=SDL2(NA1(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q3(I)*Q1(J)/D)*331.979 ET=ET+E+EQ IF(KC.EQ.2)GO TO 12 IF(D.LE.DMAX)WRITE(16,100)(AT(NA3(I),L),L=1,6), 1 (AT(NA1(J),L),L=1,6),D,E,ET 100 FORMAT(1X,6A1,'...',6A1,5X,'d =',F8.5,8X,'E =',F10.3,3X,'ET=', 1 F10.3) 12 CONTINUE DO 19 I=3,NG1-1 DO 19 J=I+1,NG1 DO 14 K=1,8 IF(NA1(J).EQ.MC(NA1(I),K))GO TO 19 DO 13 L=1,8 IF(NA1(J).EQ.MC(MC(NA1(I),K),L))GO TO 19 13 CONTINUE 14 CONTINUE IF(KC.EQ.2)GO TO 16 DO 15 L=3,NG2 IF(NA1(J).EQ.NA2(L))GO TO 19 15 CONTINUE 16 DO 17 K=1,3 X1(K)=XO1(I,K) 17 X2(K)=XO1(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 19 IF(D.GT.10.0)GO TO 19 KS(1)=SDL2(NA1(I)) KS(2)=SDL2(NA1(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q1(I)*Q1(J)/D)*331.979 ET=ET+E+EQ 19 CONTINUE IF(KC.EQ.1) GO TO 36 ELSEIF(KD.EQ.2)THEN DO 24 I=1,NG3 DO 24 J=3,NG2 DO 21 K=1,8 IF(NA2(J).EQ.MC(NA3(I),K))GO TO 24 DO 20 L=1,8 IF(NA2(J).EQ.MC(MC(NA3(I),K),L))GO TO 24 20 CONTINUE 21 CONTINUE DO 22 K=1,3 X1(K)=XO3(I,K) 22 X2(K)=XO2(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 24 IF(D.GT.10.0)GO TO 24 KS(1)=SDL2(NA3(I)) KS(2)=SDL2(NA2(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q3(I)*Q2(J)/D)*331.979 ET=ET+E+EQ 24 CONTINUE DO 30 I=3,NG1 DO 30 J=3,NG2 DO 26 K=1,8 IF(NA2(J).EQ.MC(NA1(I),K))GO TO 30 DO 25 L=1,8 IF(NA2(J).EQ.MC(MC(NA1(I),K),L))GO TO 30 25 CONTINUE 26 CONTINUE DO 27 L=3,NG2 IF(NA1(I).EQ.NA2(L))GO TO 30 27 CONTINUE DO 28 K=1,3 X1(K)=XO1(I,K) 28 X2(K)=XO2(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 30 IF(D.GT.10.0)GO TO 30 KS(1)=SDL2(NA1(I)) KS(2)=SDL2(NA2(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q1(I)*Q2(J)/D)*331.979 ET=ET+E+EQ 30 CONTINUE DO 35 I=3,NG2-1 DO 35 J=I+1,NG2 DO 32 K=1,8 IF(NA2(J).EQ.MC(NA2(I),K))GO TO 35 DO 31 L=1,8 IF(NA2(J).EQ.MC(MC(NA2(I),K),L))GO TO 35 31 CONTINUE 32 CONTINUE DO 33 K=1,3 X1(K)=XO2(I,K) 33 X2(K)=XO2(J,K) CALL DIST(X1,X2,D) IF (D.LT.0.1) GO TO 35 IF(D.GT.10.0)GO TO 35 KS(1)=SDL2(NA2(I)) KS(2)=SDL2(NA2(J)) CALL ENERPOT(KS,D,E,IND) IF(IND.EQ.1)E=0 EQ=(Q2(I)*Q2(J)/D)*331.979 ET=ET+E+EQ 35 CONTINUE 36 ENDIF RETURN END SUBROUTINE DIST(X1,X2,D) C-----This routine calculates the distance between two atoms DIMENSION X1(3),X2(3) D=0.0 DO 2 K=1,3 AA=ABS(X1(K)-X2(K)) IF(AA.GT.1E-9)GO TO 1 AA=0.0 GO TO 2 1 AA=AA*AA 2 D=D+AA D=SQRT(D) RETURN END SUBROUTINE READAT(NG,T1,T2,DT,ATT,ATR) C-----This routine reads the data defining the rotation of the fragment CHARACTER *1 ATT(4,6),MINUS(2),ATR(500,6) COMMON/INPO/IN,IO READ(IN,*)NG,T1,T2,DT READ(IN,'(4(6A1))')((ATT(I,L),L=1,6),I=1,4) DO 1 I=1,4 DO 1 K=1,3,2 MINUS(1)=ATT(I,K) MINUS(2)=ATT(I,K+1) CALL MINUMA(MINUS) ATT(I,K)=MINUS(1) 1 ATT(I,K+1)=MINUS(2) READ(IN,'(12(6A1))')((ATR(I,L),L=1,6),I=1,NG) DO 2 J=1,NG DO 2 K=1,3,2 MINUS(1)=ATR(J,K) MINUS(2)=ATR(J,K+1) CALL MINUMA(MINUS) ATR(J,K)=MINUS(1) 2 ATR(J,K+1)=MINUS(2) RETURN END SUBROUTINE ORTHOG(C,A) C-----This routine calculates the orthogonalization matrix, O(3,3) DIMENSION C(3),A(3),SN(3),CS(3) COMMON/GRP1/O(3,3) RD=0.0174532925 DO 1 I=1,3 AR=A(I)*RD SN(I)=SIN(AR) 1 CS(I)=COS(AR) CSAR=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3)) SNAR=SQRT(1.0-CSAR**2) O(1,1)=C(1) O(1,2)=C(2)*CS(3) O(1,3)=C(3)*CS(2) O(2,1)=0.0 O(2,2)=C(2)*SN(3) O(2,3)=-C(3)*SN(2)*CSAR O(3,1)=0.0 O(3,2)=0.0 O(3,3)=C(3)*SN(2)*SNAR RETURN END SUBROUTINE MINUMA(STR) C-----This routine transforms the lower case characters into upper C-----case characters, in a string of two characters CHARACTER*1 MINU(26),MAIU(26),STR(2) DATA MINU/'a','b','c','d','e','f','g','h','i','j','k','l','m', 1 'n','o','p','q','r','s','t','u','v','w','x','y','z'/ DATA MAIU/'A','B','C','D','E','F','G','H','I','J','K','L','M', 1 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/ DO 2 I=1,2 DO 1 J=1,26 IF(STR(I).NE.MINU(J)) GO TO 1 STR(I)=MAIU(J) GO TO 2 1 CONTINUE 2 CONTINUE RETURN END SUBROUTINE NORM(X) C-----This routine normalizes a 3-dimensional vector DOUBLE PRECISION X(3) S=0.0 DO 1 I=1,3 1 S=S+X(I)**2 S=SQRT(S) DO 2 I=1,3 2 X(I)=X(I)/S RETURN END SUBROUTINE NEWCOO(NG,XNO,BI,XOO,RR) C-----Modifies coordinates to have them referred to a Cartesian system C-----of axis with the origin at the first atom, with the z axis along C-----the first to the second atom, the y axis perpendicular to the z C-----axis and lying on the plane defined by the first three atoms, and C-----the x axis perpendicular to y and z. DIMENSION XNO(500,3),X1(500,3),XOO(500,3),RR(3) DOUBLE PRECISION V,XR(3),B(3,3),BI(3,3) DO 1 I=1,NG DO 1 J=1,3 1 X1(I,J)=XNO(I,J)-XNO(1,J) C-----Building up matrix B DO 2 I=1,3 XR(I)=X1(2,I) 2 RR(I)=XNO(1,I) CALL NORM(XR) DO 3 I=1,3 3 B(3,I)=XR(I) V=0 DO 4 K=1,3 4 V=V+X1(3,K)*B(3,K) DO 5 K=1,3 5 XR(K)=X1(3,K)-V*B(3,K) CALL NORM(XR) DO 6 I=1,3 6 B(2,I)=XR(I) XR(1)=B(2,2)*B(3,3)-B(3,2)*B(2,3) XR(2)=B(3,1)*B(2,3)-B(2,1)*B(3,3) XR(3)=B(2,1)*B(3,2)-B(3,1)*B(2,2) CALL NORM(XR) DO 7 I=1,3 7 B(1,I)=XR(I) DO 8 I=1,NG DO 8 J=1,3 XOO(I,J)=0.0 DO 8 K=1,3 8 XOO(I,J)=XOO(I,J)+B(J,K)*X1(I,K) CALL INVMAT3(B,BI) RETURN END SUBROUTINE ROTCOO(NG,XOO,FI,BI,XRO,RR) C-----Rotates the coordinates of a fragment: XOO --> XRO DIMENSION XOO(500,3),XRO(500,3),X1(500,3),RT(3,3),RR(3) DOUBLE PRECISION BI(3,3) C-----Building up rotation matrix RT(1,1)=COS(FI) RT(1,2)=-SIN(FI) RT(1,3)=0.0 RT(2,1)=SIN(FI) RT(2,2)=COS(FI) RT(2,3)=0.0 RT(3,1)=0.0 RT(3,2)=0.0 RT(3,3)=1.0 C-----Applying rotation matrix DO 1 I=1,NG DO 1 J=1,3 X1(I,J)=0.0 DO 1 K=1,3 1 X1(I,J)=X1(I,J)+RT(J,K)*XOO(I,K) DO 2 I=1,NG DO 2 J=1,3 XRO(I,J)=0.0 DO 2 K=1,3 2 XRO(I,J)=XRO(I,J)+BI(J,K)*X1(I,K) DO 3 I=1,NG DO 3 J=1,3 XRO(I,J)=XRO(I,J)+RR(J) 3 CONTINUE RETURN END SUBROUTINE INVMAT3(B,BI) C-----This routine inverts a 3x3 matrix DOUBLE PRECISION B(3,3),BI(3,3) BI(1,1)=B(2,2)*B(3,3)-B(2,3)*B(3,2) BI(2,1)=B(3,1)*B(2,3)-B(2,1)*B(3,3) BI(3,1)=B(2,1)*B(3,2)-B(3,1)*B(2,2) BI(1,2)=B(3,2)*B(1,3)-B(1,2)*B(3,3) BI(2,2)=B(1,1)*B(3,3)-B(3,1)*B(1,3) BI(3,2)=B(3,1)*B(1,2)-B(1,1)*B(3,2) BI(1,3)=B(1,2)*B(2,3)-B(2,2)*B(1,3) BI(2,3)=B(2,1)*B(1,3)-B(1,1)*B(2,3) BI(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1) RETURN END SUBROUTINE ENERPOT(KS,D,E,IND) C-----This routine calculates the van der Waals potential energy for a C-----non bonding contact between two atoms (only H, C, N, O, F, Si, P, C-----S, Cl, Fe, Br, Sn, I, Pb). C-----Coefficients from K.Mirsky (Computing in Crystallography, 1978) C-----for H, C, N, O, F, S, Cl; from P. Johnston et al. (J.Cryst. C-----Spectr.Res.,18, 403, 1988) for P, Fe I; from Ahmed et al. (Acta C-----Cryst. B27, 867, 1971) for Si, Sn, Pb; from Eliel et al. (Con- C-----formational Analysis, Chap.7, Interscience Publ.,1965) for Br. CHARACTER*2 ATM(14),KS(2) DIMENSION COA(14),COB(14),COC(14) DATA ATM/'H ','C ','N ','O ','F ','SI','P ','S ','CL','FE','BR', 1 'SN','I ','PB'/ DATA COA/4900,71600,42000,77700,42000,42000,525963,235000,4580, 1 70301,359377,42000,38994,42000/ DATA COB/4.29,3.68,3.78,4.18,4.15,2.95,4.0,3.49,2.262,3.18,3.484, 1 2.85,2.5,3.06/ DATA COC/29,421,259,259.4,148,1160,1531.54,2346,2980,277.8, 1 3436.3,1401,12451,912/ IF(D.EQ.0.)RETURN IND=0 DO 1 I=1,14 IF(KS(1).EQ.ATM(I)) GO TO 2 1 CONTINUE IND=1 RETURN 2 K1=I DO 3 I=1,14 IF(KS(2).EQ.ATM(I)) GO TO 4 3 CONTINUE IND=1 RETURN 4 K2=I A=COA(K1) B=COB(K1) C=COC(K1) F=(COC(K1)+COC(K2))/COC(1) IF(K1.NE.K2) GO TO 5 GO TO 6 5 A=SQRT(COA(K1)*COA(K2)) B=(COB(K1)+COB(K2))/2 C=SQRT(COC(K1)*COC(K2)) 6 E=A*EXP(-B*D)-C*D**(-6)+F*D**(-12) RETURN END SUBROUTINE TORS(R,TAU) C-----This routine calculates the torsion angles formed by the atoms C-----1-2-3-4.The angle is positive when the 1-2 bond, view down the C-----2-3 bond, will eclipse the 3-4 bond when rotates less than 180 C-----degrees in a clockwise direction ("right-hand rule",Klyne,W.& C-----Prelog,V. (1960), Experientia, 16, 521). REAL M(3,3) DIMENSION R(4,3),DLT(3,3),DST(3) DO 1 J=1,3 DLT(1,J)=R(2,J)-R(1,J) DLT(2,J)=R(2,J)-R(3,J) 1 DLT(3,J)=R(4,J)-R(3,J) DO 2 I=1,3 DST(I)=SQRT(DLT(I,1)**2+DLT(I,2)**2+DLT(I,3)**2) IF(DST(I).EQ.0.) GO TO 4 2 CONTINUE DO 3 I=1,3 DO 3 J=1,3 3 M(I,J)=DLT(I,J)/DST(I) CF1=M(1,1)*M(2,1)+M(1,2)*M(2,2)+M(1,3)*M(2,3) CF2=M(2,1)*M(3,1)+M(2,2)*M(3,2)+M(2,3)*M(3,3) IF(ABS(CF1).GT..98.OR.ABS(CF2).GT..98 )GO TO 4 SF1=SQRT(1.-CF1**2) SF2=SQRT(1.-CF2**2) P=SF1*SF2 IF(ABS(P).LT.1.E-5)GO TO 4 STAU=(M(3,1)*(M(2,2)*M(1,3)-M(1,2)*M(2,3)) 1 +M(3,2)*(M(1,1)*M(2,3)-M(1,3)*M(2,1)) 2 +M(3,3)*(M(2,1)*M(1,2)-M(1,1)*M(2,2)))/P CTAU=((M(2,2)*M(1,3)-M(1,2)*M(2,3))*( M(3,2)*M(2,3)-M(2,2)*M(3,3) 1 )+(M(2,3)*M(1,1)-M(2,1)*M(1,3))*(M(3,3)*M(2,1)-M(3,1)*M(2,3)) 2 +(M(2,1)*M(1,2)-M(1,1)*M(2,2))*(M(3,1)*M(2,2)-M(2,1)*M(3,2)))/P RR=57.295779 TAU1=ATAN2(STAU,CTAU) TAU=RR*TAU1 IF(TAU.GT.180.)TAU=TAU-360. 4 RETURN END SUBROUTINE BONDI(N) C-----This routine calculates bond distances and set up the the C-----connectivity matrix MC CHARACTER*1 ATOM(2,2),AT(500,6) CHARACTER*2 NOMAT(103),KS(2),SDL2(500) DIMENSION RG(2),IRAG(103) COMMON/GRP2/AT 2 /GRP3/MC(400,8) 3 /GRP5/XO(500,3) 4 /SMBL/NOMAT 5 /ATRD/IRAG 6 /GRG1/SDL2 7 /INPO/IN,IO DO 1 J=1,N KK=0 DO 1 K=1,N D1=0.0 DO 2 I=1,3 DR=XO(J,I)-XO(K,I) IF(ABS(DR).LT.1E-9)THEN DR=0.0 ELSE DR=DR**2 ENDIF 2 D1=D1+DR D1=SQRT(D1) IF(D1.LT.1E-3) GO TO 1 KS(1)=SDL2(J) KS(2)=SDL2(K) IF (KS(1).EQ.NOMAT(1).AND.KS(2).EQ.NOMAT(1))GO TO 1 DO 6 JJ=1,2 DO 4 I=1,103 IF(KS(JJ).EQ.NOMAT(I)) GO TO 5 4 CONTINUE WRITE(*,'(/'' There is something wrong with the atom symbols''/ 1 '' Check the atom labels''/)') GO TO 1 5 RG(JJ)=IRAG(I) 6 CONTINUE DSTC=(RG(1)+RG(2))/100. IF(D1.GT.(DSTC+0.5))GO TO 1 KK=KK+1 IF (KK.GT.8) GO TO 1 MC(J,KK)=K 1 CONTINUE RETURN END BLOCK DATA CHARACTER*2 NOMAT(103) COMMON/SMBL/NOMAT 1 /ATRD/IRAG(103) DATA NOMAT/'H ','HE','LI','BE','B ','C ','N ','O ','F ','NE', 1 'NA','MG','AL','SI','P ','S ','CL','AR','K ','CA','SC', 2 'TI','V ','CR','MN','FE','CO','NI','CU','ZN','GA','GE', 3 'AS','SE','BR','KR','RB','SR','Y ','ZR','NB','MO','TC', 4 'RU','RH','PD','AG','CD','IN','SN','SB','TE','I ','XE', 5 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB', 6 'DY','HO','ER','TM','YB','LU','HF','TA','W ','RE','OS', 7 'IR','PT','AU','HG','TL','PB','BI','PO','AT','RN','FR', 8 'RA','AC','TH','PA','U ','NP','PU','AM','CM','BK','CF', 9 'ES','FM','MD','NO','LW'/ DATA IRAG/37,93,134,90,82,77,75,73,72,131,154,130,118, 1 111,106,102,99,174,196,174,144,136,117,128,122,130,126, 2 140,138,131,126,122,119,116,114,189,211,192,162,148, 3 139,133,135,140,145,148,153,148,144,141,138,135,133, 4 209,225,198,169,165,165,164,165,166,185,161,159,159, 5 158,157,156,170,156,157,143,137,128,126,132,150,150, 6 149,148,147,146,168,150,214,230,204,188,165,150,142, 7 165,166,185,161,160,160,160,160,160,160,160/ END SUBROUTINE CONTACT(BRV,N,NR,AFR,NL,AAAA,XXO,QF) C-----Calculates the coordinates of the environment of the rotating C-----atoms CHARACTER*1 BRV,STPE(38,20),QS,AAAA(500,6),MINUS(2),LTP(7), 1 AFR(500,6),AT(500,6) DIMENSION XXO(500,3),QF(500),NS(3) COMMON/LBL1/C(3),A(3) 1 /LBL3/IC,NE,NT 2 /LBL4/STPE 3 /LBL8/XX(500,3) 4 /GRP1/O(3,3) 5 /GRP2/AT 6 /GRP5/XO(500,3) 7 /GRP10/X(500,3),Q(500) 8 /INS4/CS(3),SN(3) 9 /INPO/IN,IO DATA LTP /'P','A','B','C','I','F','R'/ PG=3.14159265 RD=57.295779513 DO 1 I=1,7 IF(BRV.EQ.LTP(I))GO TO 2 1 CONTINUE WRITE(IO,102) BRV 102 FORMAT(/' The symbol ',A1,'for lattice type is uncorrect.'/ 1 ' Symbol P has been assumed'/) BRV=LTP(1) 2 CONTINUE D3=2.0 DM=4.0 C-----Intermolecular contacts less than DM=4.0 A NL=0 CALL INTERC(N,DM,NL,BRV,NR,AFR,AAAA) C-----File of coordinates for rotener CALL ELIMRA(NL,AAAA) DO 7 I=1,NL DO 6 J=1,N DO 4 L=1,6 IF(AAAA(I,L).NE.AT(J,L))GO TO 6 4 CONTINUE DO 5 JJ=1,3 XXO(I,JJ)=0.0 DO 5 K=1,3 5 XXO(I,JJ)=XXO(I,JJ)+O(JJ,K)*XX(I,K) QF(I)=Q(J) GO TO 7 6 CONTINUE WRITE(IO,'('' There is some error in the atom labels'')') GO TO 8 7 CONTINUE 8 RETURN END SUBROUTINE INTERC(N,DM,NL,BRV,NR,AFR,AAAA) C-----This routine calculates all the interatomic contacts, C-----involving atoms of different equivalent positions in the same C-----or in adjacent cells, less than a given maximum value CHARACTER*1 STR(20),BRV,STPE(38,20),AT(500,6),ATOM(2,2), 1 AD2(6),AFR(500,6),AAAA(500,6) INTEGER LE1(24,3,3),LE(3,3) DIMENSION B(3,3),F1(24,3),F(3),FB(3) COMMON/LBL1/C(3),A(3) 1 /LBL3/IC,NE,NT 3 /LBL4/STPE 4 /LBL8/XX(500,3) 5 /INS4/CS(3),SN(3) 6 /GRP2/AT 7 /GRP10/X(500,3),Q(500) 8 /EQVP/LE1,F1,FB 9 /COOA/X1(3),X2(3) DO 1 I=1,24 DO 1 J=1,3 F1(I,J)=0.0 DO 1 K=1,3 1 LE1(I,J,K)=0 NE=NE+1 IF(NE.EQ.1)GO TO 4 DO 3 I=2,NE DO 2 L=1,20 2 STR(L)=STPE(I-1,L) CALL DECEQP(STR,LE,F) DO 3 J=1,3 F1(I,J)=F(J) LE1(1,J,J)=1 DO 3 K=1,3 3 LE1(I,J,K)=LE(J,K) GO TO 6 4 DO 5 J=1,3 5 LE1(1,J,J)=1 6 CONTINUE DO 7 I=1,3 DO 7 J=1,3 7 B(I,J)=0.0 KB=0 IF(BRV.EQ.'P') GO TO 13 IF(BRV.NE.'A') GO TO 8 B(1,2)=1./2. B(1,3)=1./2. KB=1 GO TO 13 8 IF(BRV.NE.'B') GO TO 9 B(1,1)=1./2. B(1,3)=1./2. KB=1 GO TO 13 9 IF(BRV.NE.'C') GO TO 10 B(1,1)=1./2. B(1,2)=1./2. KB=1 GO TO 13 10 IF(BRV.NE.'I') GO TO 11 B(1,1)=1./2. B(1,2)=1./2. B(1,3)=1./2. KB=1 GO TO 13 11 IF(BRV.NE.'R') GO TO 12 B(1,1)=1./3. B(2,2)=1./3. B(2,3)=1./3. B(1,2)=2*B(1,1) B(1,3)=2*B(1,1) B(2,1)=2*B(1,1) KB=2 GO TO 13 12 B(1,2)=1./2. B(1,3)=1./2. B(2,1)=1./2. B(2,3)=1./2. B(3,1)=1./2. B(3,2)=1./2. KB=3 13 CONTINUE DO 22 II=1,NR DO 21 I=1,N DO 14 L=1,6 IF(AT(I,L).NE.AFR(II,L))GO TO 21 14 CONTINUE DO 20 J=1,N DO 15 L=1,6 15 AD2(L)=AT(J,L) DO 16 K=1,3 X1(K)=X(I,K) 16 X2(K)=X(J,K) DO 17 L=1,3 17 FB(L)=0.0 CALL DISCON(AD2,DM,NL,AAAA) IF(KB.EQ.0)GO TO 20 DO 19 L=1,KB DO 18 K=1,3 18 FB(K)=B(L,K) CALL DISCON(AD2,DM,NL,AAAA) 19 CONTINUE 20 CONTINUE 21 CONTINUE 22 CONTINUE RETURN END SUBROUTINE DISCON(AD2,DM,NL,AAAA) C-----Set data and calculates contacts involving atoms translated in C-----equivalent positions CHARACTER*1 STPE(38,20),AD2(6),AAAA(500,6) INTEGER LE1(24,3,3),FT(3),LE(3,3) DIMENSION F1(24,3),FB(3),F(3) COMMON/COOA/X1(3),X2(3) 1 /LBL3/IC,NE,NT 2 /LBL1/C(3),A(3) 3 /INS4/CS(3),SN(3) 4 /EQVP/LE1,F1,FB 5 /LBL4/STPE 6 /LBL8/XX(500,3) 7 /INPO/IN,IO KM=2*NT+1 DO 3 I=1,NE DO 3 I1=1,KM FT(1)=I1-NT-1 DO 3 I2=1,KM FT(2)=I2-NT-1 DO 3 I3=1,KM FT(3)=I3-NT-1 DO 1 J=1,3 F(J)=F1(I,J)+FB(J)+FT(J) DO 1 K=1,3 1 LE(J,K)=LE1(I,J,K) CALL DSTCHB(DM,AD2,LE,F,NL,AAAA) IF(IC.NE.-1) GO TO 3 DO 2 J=1,3 F(J)=-F1(I,J)+FB(J)+FT(J) DO 2 K=1,3 2 LE(J,K)=-LE1(I,J,K) CALL DSTCHB(DM,AD2,LE,F,NL,AAAA) 3 CONTINUE RETURN END SUBROUTINE DSTCHB(DM,AD2,LE,F,NL,AAAA) C-----This routine calculates interatomic distances involving C-----atoms translated in equivalent positions CHARACTER*1 AD2(6),AAAA(500,6) INTEGER LE(3,3) DIMENSION F(3),X(3),D(3) COMMON/COOA/X1(3),X2(3) 1 /LBL3/IC,NE,NT 2 /LBL1/C(3),A(3) 3 /INS4/CS(3),SN(3) 4 /LBL8/XX(500,3) 5 /INPO/IN,IO DO 1 I=1,3 X(I)=F(I) DO 1 J=1,3 1 X(I)=X(I)+LE(I,J)*X2(J) DO 2 I=1,3 2 IF(ABS(X(I)-X2(I)).GT.1E-4)GO TO 3 GO TO 8 3 DA=0 DO 4 I=1,3 D(I)=X1(I)-X(I) D(I)=D(I)*C(I) 4 DA=DA+D(I)**2 IF(DA.LT.1E-2)GO TO 8 DA=DA+2*(D(1)*D(2)*CS(3)+D(1)*D(3)*CS(2)+D(2)*D(3)*CS(1)) DA=SQRT(DA) IF(DA.GT.DM) GO TO 8 C-----Check for duplicates DO 11 I=1,NL DO 9 L=1,6 IF(AD2(L).NE.AAAA(NL,L))GO TO 11 9 CONTINUE DO 10 K=1,3 IF(X(K).NE.XX(NL,K))GO TO 11 10 CONTINUE GO TO 8 11 CONTINUE NL=NL+1 IF(NL.LE.500)GO TO 5 WRITE(IO,101) 101 FORMAT(' There are too many atoms whose environment is ', 1 'requested'/' Please, reduce their number'/) NL=500 GO TO 8 5 DO 6 K=1,3 6 XX(NL,K)=X(K) DO 7 L=1,6 7 AAAA(NL,L)=AD2(L) 8 RETURN END SUBROUTINE DECEQP(STR,LE,F) C-----This routine decodes the equivalent positions CHARACTER*1 STR(20),COO(13) INTEGER LE(3,3) DIMENSION F(3) COMMON/INPO/IN,IO DATA COO/'X','Y','Z',',','/','+','-','1','2','3','4','5','6'/ DO 1 I=1,3 F(I)=0.0 DO 1 J=1,3 1 LE(I,J)=0 N=0 DO 3 I=1,20 IF(STR(I).EQ.' '.AND.STR(I+1).EQ.' ')GO TO 4 IF(STR(I).NE.COO(4))GO TO 3 N=N+1 IF(N.EQ.1)NA=I IF(N.EQ.2)NB=I 3 CONTINUE 4 DO 17 K=1,3 GO TO (5,6,7),K 5 N1=1 N2=NA-1 GO TO 8 6 N1=NA+1 N2=NB-1 GO TO 8 7 N1=NB+1 N2=20 8 KS=1 KD=0 DO 16 I=N1,N2 IF(STR(I).EQ.' '.AND.STR(I+1).EQ.' ')GO TO 18 IF(STR(I).EQ.' ')GO TO 16 DO 9 J=1,13 IF(STR(I).EQ.COO(J))GO TO 10 9 CONTINUE GO TO 16 10 IF(J.EQ.4)GO TO 16 IF(J.EQ.7)KS=-1 IF(J.EQ.5)KD=1 IF(J.GT.3)GO TO 12 DO 11 L=1,3 11 IF(K.EQ.L)LE(L,J)=1*KS KS=1 GO TO 16 12 IF(J.LT.8)GO TO 16 IF(KD.EQ.1)GO TO 14 DO 13 L=1,3 13 IF(K.EQ.L)F(L)=(J-7)*KS KS=1 GO TO 16 14 DO 15 L=1,3 IF(K.EQ.L)F(L)=F(L)/(J-7) 15 CONTINUE KD=0 16 CONTINUE 17 CONTINUE 18 RETURN END SUBROUTINE ELIMRA(NL,AAAA) C-----Eliminate repeated atoms CHARACTER*1 AAAA(500,6),AAA1(500,6) DIMENSION X1(500,3) COMMON/LBL8/XX(500,3) DO 6 I=1,NL-1 DO 5 J=I+1,NL DO 1 L=1,6 IF(AAAA(I,L).NE.AAAA(J,L))GO TO 5 1 CONTINUE DO 2 K=1,3 IF(ABS(XX(I,K)-XX(J,K)).GT.0.0001)GO TO 5 2 CONTINUE DO 3 L=1,6 3 AAAA(J,L)=' ' DO 4 K=1,3 4 XX(J,K)=0.0 5 CONTINUE 6 CONTINUE NN=0 DO 11 I=1,NL DO 7 L=1,6 IF(AAAA(I,1).NE.' ') GO TO 8 7 CONTINUE GO TO 11 8 NN=NN+1 DO 9 L=1,6 9 AAA1(NN,L)=AAAA(I,L) DO 10 K=1,3 10 X1(NN,K)=XX(I,K) 11 CONTINUE DO 13 I=1,NN DO 12 L=1,6 12 AAAA(I,L)=AAA1(I,L) DO 13 K=1,3 13 XX(I,K)=X1(I,K) NL=NN RETURN END SUBROUTINE ADD(N,NL,AAAA,XO3,XXO,Q3,QF,NG3) CHARACTER*1 AT(500,6),AAAA(500,6) DIMENSION XO3(500,3),XXO(500,3),Q3(500),QF(500) COMMON/GRP2/AT 1 /GRP6/NA1(500),NA2(500),NA3(500) 2 /INPO/IN,IO DO 2 I=1,NL DO 1 K=1,3 1 XO3(I+NG3,K)=XXO(I,K) Q3(I+NG3)=QF(I) NA3(I+NG3)=I+N DO 2 L=1,6 2 AT(I+N,L)=AAAA(I,L) NG3=NG3+NL N=N+NL RETURN END