PROGRAM CART C ************************************************************************** C ************************************************************************** C THIS PROGRAM USES INPUT FROM CARTIN.DAT. ITS OUTPUT IS TO CARTOUT.DAT. C THIS IS A MAJOR MODIFICATION OF PICKETT'S PROGRAM CART IN WHICH MOST OF C THE JUNK THAT DOES NOT APPLY TO OUR WORK HAS BEEN REMOVED AND OPTION C CODING HAS BEEN SIMPLIFIED. FOR DOCUMENTATRION SEE SORENSON'S (7-87). C THE ATOMIC WEIGHT TABLE WILL NOT BE PRINTED OUT. WHEN DOING INTERNAL C ROTATION, B+C AND RAY'S PARAMETER WILL BE PRINTED OUT. IF IOPT(4) IS C SET TO 9, CERTAIN DATA WILL BE LEFT OUT IN ALL ITERATIONS AFTER THE C FIRST AND B+C WILL NOT BE PRINTED FOR KAPPA LESS NEGATIVE THAN -0.60. C ************************************************************************** C ************************************************************************** IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) INTEGER ATNO REAL*4 SYMB LOGICAL CALC,TAPE,NOTAPE,NODER,NOISO,RESET DIMENSION ATNO(20),MSNO(20),SYMB(20),TITLE(9),DA(3),DB(3),DC(3), A IQ(3),SCL(3) COMMON/POT/IAF(50),JAF(50),NQ,NNATM COMMON XYZ(50,3),Q(50,3),F(50,3),WT(20),NATM,NO(51),MASS(50), B NA(50),NB(50),NC(50),NRDN,NZ,IOPT(4) C ************************************************************************ C ************************************************************************ C IF YOU WISH TO ADD TO OR CORRECT THE ATOMIC WEIGHT TABLE YOU WILL HAVE TO C ALTER LINES 12, 14, 15, 40 AND 46 (CART NUMBERS) AND MAKE PROPER CHANGES C IN THE TABLE (LINES 17 THROUGH 36). LOTSA LUCK! C ************************************************************************ C ************************************************************************ DATA SYMB/'C 12','H 1','O 16','N 14','F 19','H 2','O 18','C 13' D ,'N 15','CL35','CL37','BR79','BR81','s 32','s 34','p 31','I ', D 'DUMY','DUMY','DUMY'/ DATA ATNO/6,1,8,7,9,1,8,6,7,17,17,35,35,16,16,15,53,3*0/ DATA MSNO/12,1,16,14,19,2,18,13,15,35,37,79,81,32,34,31,127,3*0/ DATA SCL,DR/0.001,3*0.0174532925199D 00/ OPEN (UNIT=5, FILE= 'CARTIN.DAT',STATUS= 'OLD') OPEN (UNIT=6, FILE= 'CARTOUT.DAT',STATUS= 'NEW') WT( 1)=12. WT( 2)= 1.00782 WT( 3)=15.99491 WT( 4)=14.00307 WT( 5)=18.99840 WT( 6)= 2.01409 WT( 7)=17.99915 WT( 8)=13.00335 WT( 9)=15.00011 WT(10)=34.96885 WT(11)=36.96582 WT(12)=78.9183 WT(13)= 80.9163 WT(14)= 31.9722 WT(15)= 33.9688 WT(16)= 30.9744 WT(17)= 126.9004 WT(18)= 0. WT(19)= 0. WT(20)= 0. NO(51)=0 NATM=0 N=17 C ************************************************************************** C ************************************************************************** READ(5,100) J,TITLE 100 FORMAT(I5,9A8) IF(J.LE.0) GO TO 1 C EXTEND ATOMIC WEIGHT LIBRARY N=N+J READ(5,101)(SYMB(I),ATNO(I),MSNO(I),WT(I),I=18,N) 101 FORMAT(A4,1X,2I5,F20.5) C EACH STEP STARTS HERE 1 READ(5,1001) NZ,TITLE 1001 FORMAT(I5,9A8) IF(NZ) 3,2,4 2 CALL EXIT C ************************************************************************** C PROCEEDURES FOR STEP IN WHICH DATA IS USED FROM LAST STEP 3 IF(NATM.LE.0) GO TO 2 N=-NZ READ (5,103) IFDER,M ,IOPT 103 FORMAT(6I5) NOISO=.FALSE. DO 5 I=1,N READ(5,104) M,IW,JW,IA,IB,IC,AT,BT,CT 104 FORMAT(6I3,2X,6F10.3) DO 6 K=1,NATM IF(NO(K).NE.M) GO TO 6 IF(IA.EQ.NO(NA(K)).AND.IB.EQ.NO(NB(K)).AND.IC.EQ.NO(NC(K)))GO TO 7 6 CONTINUE if(iopt(4).eq.9)go to 1004 9 WRITE(6,105)M,IW,JW,IA,IB,IC,AT,BT,CT 1004 continue 105 FORMAT(32H1 INPUT CARD ERROR -- CARD WAS 6I3,2X,6F10.4) GO TO 2 7 IF(DABS(Q(K,1)-AT).LT.0.000001) GO TO 70 Q(K,1)=AT NOISO=.TRUE. 70 IF(DABS(Q(K,2)/DR-BT).LT.0.000001) GO TO 71 Q(K,2)=BT*DR NOISO=.TRUE. 71 IF(DABS(Q(K,3)/DR-CT).LT.0.000001) GO TO 72 Q(K,3)=CT*DR NOISO=.TRUE. 72 IF(IW.EQ.ATNO(MASS(K)).AND.JW.EQ.MSNO(MASS(K))) GO TO 5 DO 8 J=1,20 IF(IW.NE.ATNO(J)) GO TO 8 IF(JW.NE.MSNO(J)) GO TO 8 MASS(K)=J GO TO 5 8 CONTINUE 28 WRITE(6,112) IW,JW 112 FORMAT(35H1AN ATOMIC WEIGHT FOR ATOMIC NUMBER I4,14H AND MASS NUMB 1 2HER,I4, 22H IS NOT IN THE LIBRARY ) GO TO 2 5 CONTINUE IF(NOISO) GO TO 11 C SPECIAL ACTION WHEN ONLY MASSES ARE CHANGED CALC=(NODER.AND.IFDER.NE.0).OR.NOTAPE RESET=.NOT.CALC TAPE=.FALSE. GO TO 10 C ACTION WHEN COORDINATES ARE CHANGED FROM LAST STEP 11 CALC=.TRUE. NODER=IFDER.EQ.0 TAPE=IFSAVE.NE.0 NOTAPE=.NOT.TAPE RESET=TAPE 10 DO 13 I=1,NNATM Q(I,2)=Q(I,2)/DR 13 Q(I,3)=Q(I,3)/DR GO TO 12 C ************************************************************************** C PROCEEDURES WHEN STEP IS COMPLETELY NEW 4 READ (5,103) IFDER,NRDN,IOPT CALC=.TRUE. NODER=IFDER.EQ.0 RESET=TAPE NATM=NZ NNATM=NZ+NRDN NQ=0 DO 14 I=1,NNATM READ(5,104) NO(I),IW,JW,IA,IB,IC,(Q(I,J),J=1,6) DO 15 J=1,20 IF(IW.NE.ATNO(J)) GO TO 15 IF(JW.NE.MSNO(J)) GO TO 15 MASS(I)=J GO TO 16 15 CONTINUE GO TO 28 16 J=51 IF(I.EQ.1) GO TO 46 DO 45 J=1,I IF(NO(J).EQ.IA) GO TO 46 45 CONTINUE 29 WRITE(6,105)NO(I),IW,JW,IA,IB,IC,(Q(I,J),J=1,6) GO TO 2 46 NA(I)=J J=51 IF(I.LE.2) GO TO 48 DO 47 J=1,I IF(NO(J).EQ.IB) GO TO 48 47 CONTINUE GO TO 29 48 NB(I)=J J=51 IF(I.LE.3) GO TO 43 DO 49 J=1,I IF(NO(J).EQ.IC) GO TO 43 49 CONTINUE GO TO 29 43 NC(I)=J DO 17 J=1,3 IF(F(I,J)) 18,17,19 18 F(I,J)=10000. 19 NQ=NQ+1 IAF(NQ)=I JAF(NQ)=J F(I,J)=1./DSQRT(F(I,J)) 17 CONTINUE 14 CONTINUE 12 continue if(iopt(4).eq.9)write(6,1006)title 1006 format(1h1,9a8) if(iopt(4).eq.0)write(6,106)title,natm,nrdn 106 FORMAT(1H1,9A8///19H THIS MOLECULE HAS ,I2,19H ATOMS. THERE ARE EI2,30H EXTRA DEFINITIONS OF POSITION ,///,18H INPUT COORDINATES) DO 20 I=1,NNATM J=MASS(I) IQ(1)=NO(NA(I)) IQ(2)=NO(NB(I)) IQ(3)=NO(NC(I)) IF(IOPT(4).EQ.9)GO TO 1007 WRITE(6,107) NO(I),SYMB(J),IQ(1),IQ(2),IQ(3),(Q(I,J),J=1,6) 1007 CONTINUE 107 FORMAT(I10,2X,A4,5X,3I5,6F12.4) Q(I,2)=Q(I,2)*DR 20 Q(I,3)=Q(I,3)*DR IF(TAPE) WRITE(6,108) 108 FORMAT(50H TAPE UNIT 1 WILL BE USED FOR OUTPUT IN THIS STEP ) CONSTRUCT CALL ABCXYZ(1,CALC) A=0. B=0. C=0. CALL MOMENT(A,B,C,1) QQ=0. I=0 WRITE(7,109) I,A,B,C,QQ,TITLE(1) 109 FORMAT(I5,4F15.6,7X,A8) DO 26 I=1,NATM QQ=WT(MASS(I)) WRITE(7,109)NO(I),(XYZ(I,J),J=1,3),QQ,TITLE(1) 26 CONTINUE CALCULATE OPTIONS IF (IOPT(1).NE.0)CALL DIPOLE(1) IF(IOPT(3).GT.0) CALL ROTOR(A,B,C) IF(IFDER.LE.0) GO TO 30 CALCULATE IF(IOPT(4).EQ.9)GO TO 1002 WRITE(6,110) 1002 CONTINUE 110 FORMAT(////,50H DERIVATIVES OF ROTATIONAL CONSTANTS WITH RESPECT 1 20HTO INPUT COORDINATES ,/,35H UNITS ARE MEGAHERTZ PER DEGREE 120H OR MILLI-ANGSTROM ,/,22H ATOM NO. DISTANCE A,10X,1HB,10X, 2 1HC,6X,9HTHETA A,10X,1HB,10X,1HC,8X,7HPHI A,10X,1HB,10X,1HC) DO 21 I=2,NNATM DO 22 J=1,3 DA(J)=0. DB(J)=0. DC(J)=0. IF(J.GE.I) GO TO 22 QQ=Q(I,J) Q(I,J)=QQ+SCL(J) CALL ABCXYZ(-1,CALC) Q(I,J)=QQ CALL MOMENT(A,B,C,-1) DA(J)=A DB(J)=B DC(J)=C CALCULATE DERIVATIVE OPTIONS IF (IOPT(2).GT.0)CALL DIPOLE(0) 22 CONTINUE IQ(1)=NO(NA(I)) IQ(2)=NO(NB(I)) IQ(3)=NO(NC(I)) WRITE(6,111)NO(I),(IQ(K),DA(K),DB(K),DC(K),K=1,3) 111 FORMAT(I8,2X,3(I4,3F11.3)) 21 CONTINUE IF (IOPT(2).GT.0)CALL DIPOLE(-1) 30 IF(RESET) REWIND 1 GO TO 1 END SUBROUTINE COORDS IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N) DIMENSION VECTR(2,3),TRANS(3,3),PRCOR(3) COMMON COORD(50,3),DIST(50),THETA(50),PHI(50),F(150),WT(20),NNATM, C NAR(51),MASS(50),NA(50),NB(50),NC(50),NR,NZ,IOPT(4) C THE COORDINATES OF ATOM 1 ARE SET AS THE ORIGIN COORD(1,1)=0.0 COORD(1,2)=0.0 COORD(1,3)=0.0 C THE X COORDINATE OF ATOM 2 IS CALCULATED COORD(2,1)=DIST(2) COORD(2,2)=0.0 COORD(2,3)=0.0 C THE X AND Y COORDINATE OF ATOM 3 ARE CALCULATED IA=NA(3) IB=0 IC=0 CTH=DCOS(THETA(3)) IF(IA.EQ.2) CTH=-CTH COORD(3,1)=COORD(IA,1) +DIST(3)*CTH COORD(3,2)=COORD(IA,2)+DIST(3)*DSQRT(1.-CTH*CTH) COORD(3,3)=0.0 C THE COORDINATES OF ALL THE OTHER ATOMS ARE CALCULATED NOATM=NNATM+NR IF(NOATM-4)23,9,9 9 DO 19 ILOOP=4,NOATM IAP=IA IBP=IB ICP=IC IA=NA(ILOOP) IB=NB(ILOOP) IC=NC(ILOOP) C THE TRANSFORMATION MATRIX IS OBTAINED IF(IA-IAP)12,10,12 10 IF(IB-IBP)12,11,12 11 IF(IC-ICP)12,20,12 C THE COMPONENTS OF THE X PRIME DIRECTION ARE OBTAINED 12 DO 13 I=1,3 VECTR(1,I)=-COORD(IA,I)+COORD(IB,I) 13 VECTR(2,I)=-COORD(IA,I)+COORD(IC,I) RIJ=DSQRT(VECTR(1,1)**2+VECTR(1,2)**2+VECTR(1,3)**2) RIL=0. SCALE=0.0 DO 14 I=1,3 TRANS(I,1)=VECTR(1,I)/RIJ RIL=RIL+VECTR(2,I)**2 14 SCALE=SCALE+VECTR(2,I)*TRANS(I,1) C THE COMPONENTS OF THE Y PRIME DIRECTION ARE OBTAINED RIK=RIL-SCALE**2 IF(RIK.GT.1.0D-08 ) GO TO 15 DO 17 I=1,3 17 COORD(ILOOP,I)=COORD(IA,I)-DIST(ILOOP)*TRANS(I,1) IA=0 IB=0 IC=0 GO TO 19 15 RIK=DSQRT(RIK) DO 16 I=1,3 16 TRANS(I,2)=(VECTR(2,I)-SCALE*TRANS(I,1))/RIK C THE COMPONENTS OF THE Z DIRECTION ARE OBTAINED TRANS(1,3)=TRANS(2,1)*TRANS(3,2)-TRANS(3,1)*TRANS(2,2) TRANS(2,3)=TRANS(3,1)*TRANS(1,2)-TRANS(1,1)*TRANS(3,2) TRANS(3,3)=TRANS(1,1)*TRANS(2,2)-TRANS(2,1)*TRANS(1,2) C THE SPHERICAL COORDINATES ARE TRANSFORMED TO CARTESIAN COORDINATES C IN THE PRIME SYSTEM 20 IF(ILOOP-NNATM)24,24,25 24 IF(DIST(ILOOP)) 7,8,8 7 TMP=DIST(ILOOP)*DCOS(THETA(ILOOP)) CTH=4.*RIL CTH=SCALE/DSQRT(CTH) STH=TMP*DSQRT(.5-CTH) CTH=TMP*DSQRT(.5+CTH) SB=DSIN(PHI(ILOOP)) CB=DCOS(PHI(ILOOP)) PRCOR(1)=CTH*CB-STH*SB PRCOR(2)=STH*CB+CTH*SB PRCOR(3)=-DIST(ILOOP)*DSIN(THETA(ILOOP)) GO TO 6 8 CTH=DCOS(THETA(ILOOP)) TMP=DIST(ILOOP)*DSQRT(1.-CTH*CTH) PRCOR(1)=DIST(ILOOP)*CTH PRCOR(2)=TMP*DCOS(PHI(ILOOP)) PRCOR(3)=TMP*DSIN(PHI(ILOOP)) C THE PRIME COORDINATES ARE TRANSLATED AND ROTATED TO THE COORDINATE C SYSTEM DESIRED 6 DO 22 I=1,3 SUM=0.0 DO 21 J=1,3 21 SUM=SUM+TRANS(I,J)*PRCOR(J) 22 COORD(ILOOP,I)=COORD(IA,I)+SUM 19 CONTINUE 23 RETURN C SPECIAL REVERSE CALCULATION OF DIST THETA AND PHI 25 DO 30 I=1,ILOOP IF(NAR(ILOOP).EQ.NAR(I)) GO TO 31 30 CONTINUE 31 II=I RIJ=0. DO 26 I=1,3 SUM=0. DO 27 J=1,3 27 SUM=SUM+TRANS(J,I)*(COORD(II,J)-COORD(IA,J)) RIJ=RIJ+SUM*SUM 26 PRCOR(I)=SUM RIJ=DSQRT(RIJ) IF(DIST(ILOOP)) 28,29,29 28 STH=PRCOR(3)/RIJ COORD(ILOOP,1)=-RIJ COORD(ILOOP,2)=ASIND(STH) IF(DABS(STH) .GT. 0.99999999999999D 00) GO TO 32 CTH=4.*RIL STH=-DSQRT(0.5-CTH) CTH=-DSQRT(0.5+CTH) SB=PRCOR(2)*CTH-PRCOR(1)*STH CB=PRCOR(1)*CTH+PRCOR(2)*STH COORD(ILOOP,3)=DATAN2(SB,CB) GO TO 19 29 COORD(ILOOP,1)=RIJ STH=PRCOR(1)/RIJ COORD(ILOOP,2)=ACOSD(STH) IF(DABS(STH) .GT. 0.999999999999D 00) GO TO 32 COORD(ILOOP,3)=DATAN2(PRCOR(3),PRCOR(2)) GO TO 19 32 COORD(ILOOP,3)=PHI(ILOOP) GO TO 19 END SUBROUTINE ABCXYZ(IFLG,CALC) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) INTEGER*4 ITAPE2,IFLG LOGICAL CALC COMMON X(50,3),Q(50,3),F(50,3),W(20),NOATM,NAR(51),MASS(50),NA(50) C ,NB(50),NC(50),NRDN,NZ,IOPT(4) COMMON/POT/ IFA(50),JFA(50),NR,NNATM DIMENSION DEL(50),DPL(50),Z(50),ZZ(50),FF(50),A(25,25) COMMON/DUMCOM/A,AXZD(1375) if(iopt(4).eq.9)itape2=7 if(iopt(4).eq.0)itape2=6 IF(CALC) GO TO 100 READ(1)((X(I,J),J=1,3),I=1,NOATM) RETURN CHECK FOR REDUNDANT ATOMS 100 IF(NRDN) 1,1,2 1 CALL COORDS IF(IFLG) 24,24,25 2 IF(IFLG ) 4,4,17 4 DO 6 I=1,NR DEL(I)=0. Z(I)=Q(IFA(I),JFA(I)) 6 Q(IFA(I),JFA(I))=DPL(I) ERM=ERMX VZL=VZP GO TO 10 17 NQ=0 ERMX=0. DO 3 I=1,NR IF(IFA(I).LE.NOATM) NQ=I DEL(I)=0. FF(I)=F(IFA(I),JFA(I)) ERMX=ERMX+(1.E-05/FF(I))**2 3 Z(I)=Q(IFA(I),JFA(I)) ERM=ERMX VZL=0. NQ1=NQ+1 RL=2. STEPP=1. ITR=0 IF(NQ.GT.0.AND.NQ.LT.NR) GO TO 10 NRDN=0 GO TO 1 CALCULATE PRESENT CONSTRAINT POTENTIAL ENERGY 10 CALL COORDS VZ=0. DO 11 I=1,NQ TB=(Z(I)-Q(IFA(I),JFA(I)))/FF(I) ZZ(I)=TB 11 VZ=VZ+TB*TB DO 15 I=NQ1,NR IA=IFA(I) IB=JFA(I) TB=Z(I)-X(IA,IB) IF(IB .EQ. 1 .OR. DABS(TB) .LT. 3.1415926D 00) GO TO 150 TMP = 6.2831853071796D 00 X(IA,IB)=X(IA,IB)+DSIGN(TMP,TB) TB=Z(I)-X(IA,IB) 150 TB=TB/FF(I) ZZ(I)=TB 15 VZ=VZ+TB*TB VZ=0.5*VZ IF(DABS(VZ-VZL).LT.ERM) GO TO 20 CALCULATE DERIVATIVES OF DEPENDENT COORDINATES DO 120 K=1,NQ IA=IFA(K) IB=JFA(K) QQ=Q(IA,IB) Q(IA,IB)=QQ-0.0001 CALL COORDS Q(IA,IB)=QQ+0.0001 IJ=0 DO 13 I=NQ1,NR IJ=IJ+1 13 A(IJ,K)=X(IFA(I),JFA(I)) CALL COORDS IJ=0 DO 12 I=NQ1,NR IJ=IJ+1 12 A(IJ,K)=5000.*FF(K)*(X(IFA(I),JFA(I))-A(IJ,K))/FF(I) 120 Q(IFA(K),JFA(K))=QQ CALCULATE CHANGES IN COORDINATES CALL SOLVE (ZZ,NQ,NR) CS=0. R=0. ERR=0. DO 14 I=1,NQ TA=ZZ(I)*FF(I) CS=TA*DEL(I)+CS DEL(I)=TA TA=DABS(TA) IF(TA.GT.ERR) ERR=TA 14 R=R+TA*TA R=DSQRT(R) STEP=1. IF(ERR.GT.0.00006.AND.IFLG.GT.0) 1STEP=STEPP+STEPP*CS/(RL*(RL+R+R)) RL=R STEPP=STEP DO 16 I=1,NQ IA=IFA(I) IB=JFA(I) 16 Q(IA,IB)=Q(IA,IB)+STEP*DEL(I) IF(IFLG.LT.0) GO TO 1 VZL=VZ ITR=ITR+1 WRITE(ITAPE2,961) ITR,STEP,VZ,(IFA(I),JFA(I),DEL(I),I=1,NQ) 961 FORMAT(/, 10H ITERATION ,I4,12H STEP SIZE = E15.5,5H VZ =,D15.5,/, 1 (6(2I5,F10.5))) IF(ITR-50) 10,20,20 20 IF(IFLG ) 24,24,18 COMPLETE PREPARATION FOR DERIVATIVES AND DO PRINT OUT 18 WRITE(ITAPE2,90) 90 FORMAT(//22H STRUCTURE COORDINATES) DO 23 I=1,NNATM M=NAR(I) J=NAR(NA(I)) K=NAR(NB(I)) L=NAR(NC(I)) MM=MASS(I) IF(I.LE.NOATM) GO TO 22 R=X(I,1) ERM=X(I,2)*57.29578 ERR=X(I,3)*57.29578 GO TO 23 22 R=Q(I,1) ERM=Q(I,2)*57.29578 ERR=Q(I,3)*57.29578 23 WRITE(ITAPE2,91) M,J,K,L,R,ERM,ERR,W(MM) 91 FORMAT(4I5,F10.5,2F9.2,F13.7) DO 21 I=1,NR 21 DPL(I)=Q(IFA(I),JFA(I)) VZP=VZ 25 WRITE(ITAPE2,960)(NAR(I),(X(I,J),J=1,3),I=1,NOATM) 960 FORMAT(////12X,95HTHE COORDINATES FOR THE ATOMS IN A FRAME DEFINED 1 BY ATOMS 1,2,AND 3 AS ATOMS A,B,C RESPECTIVELY///8X,8HATOM NO.,19 2X,1HX,23X,1HY,23X,1HZ,//(9X,I5,13X,F18.7,6X,F18.7,6X,F18.7)) 24 IF(NRDN) 30,30,29 29 DO 28 I=1,NR 28 Q(IFA(I),JFA(I))=Z(I) 30 RETURN END SUBROUTINE SOLVE ( B,N,M) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) DIMENSION AA(25,25),B(50),A(50,26) COMMON/DUMCOM/AA,A,DUM(75) NN=N+1 DO 10 I=1,M 10 A(I,NN)=B(I) DO 1 I=NN,M II=I-N DO 1 J=1,N 1 A(I,J)=AA(II,J) I=1 DO 2 II=2,NN DO 3 J=I,N 3 A(I,J)=0. R=1. DO 40 J=NN,M S=A(J,I) IF(DABS(S).LT.1.E-15*R) GO TO 40 C=R R=DSQRT(R*R+S*S) S=S/R C=C/R DO 4 K=II,NN T=A(J,K) A(J,K)=T*C-A(I,K)*S 4 A(I,K)=T*S+A(I,K)*C 40 CONTINUE DO 5 K=II,NN 5 A(I,K)=A(I,K)/R B(I)=A(I,NN) 2 I=II II=N I=N 6 IF(I.EQ.1) RETURN I=I-1 DO 7 J=II,N 7 B(I)=B(I)-A(I,J)*B(J) II=I GO TO 6 END SUBROUTINE MOMENT(A,B,C,IFLG) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) INTEGER*4 ITAPE2,IFLG DIMENSION COORDX(50),COORDY(50),COORDZ(50),CMASS(3),T(3,3),COSXY( 1 3,3),CMCOR(50,3),PACOR(50,3),ROT(3) COMMON COORD(50,3),Q(300),WTATOM(20),NOATM,NO(51),MASS(50),NA(153) C,IOPT(4) EQUIVALENCE (COORDX,COORD(1)),(COORDY,COORD( 51)),(COORDZ,COORD(10 11)),(COORD,CMCOR),(COORD,PACOR) ITAPE2=6 IDER=0 ITOP=0 ITOPA=0 25 WTMOL=0.0 T(1,1)=0.0 T(2,2)=0.0 T(3,3)=0.0 T(1,2)=0.0 T(1,3)=0.0 T(2,3)=0.0 CMASSX=0.0 CMASSY=0.0 CMASSZ=0.0 C CALCULATE CENTER OF MASSES AND DO 50 I=1,NOATM M=MASS(I) TEMPX1= WTATOM(M)*COORDX(I) TEMPY1= WTATOM(M)*COORDY(I) TEMPZ1= WTATOM(M)*COORDZ(I) CMASSX= CMASSX + TEMPX1 CMASSY= CMASSY + TEMPY1 CMASSZ= CMASSZ + TEMPZ1 TEMPX2= TEMPX1*COORDX(I) TEMPY2= TEMPY1*COORDY(I) TEMPZ2= TEMPZ1*COORDZ(I) T(1,1)= T(1,1)+ TEMPY2 + TEMPZ2 T(2,2)= T(2,2)+ TEMPX2+ TEMPZ2 T(3,3)= T(3,3)+ TEMPX2 + TEMPY2 T(1,2)= T(1,2)- TEMPX1 * COORDY(I) T(1,3)= T(1,3) - TEMPX1 * COORDZ(I) T(2,3)= T(2,3) - TEMPY1 * COORDZ(I) 50 WTMOL= WTMOL + WTATOM(M) CMASSX = CMASSX/WTMOL CMASSY = CMASSY/WTMOL CMASSZ = CMASSZ/WTMOL CMASX2= CMASSX * CMASSX CMASY2= CMASSY * CMASSY CMASZ2 = CMASSZ*CMASSZ T(1,1)= T(1,1)- WTMOL * (CMASY2 + CMASZ2) T(2,2)= T(2,2)- WTMOL * (CMASX2 + CMASZ2) T(3,3)= T(3,3)- WTMOL * (CMASX2 + CMASY2) T(1,2)= T(1,2) + WTMOL * CMASSX * CMASSY T(1,3)= T(1,3) + WTMOL * CMASSX* CMASSZ T(2,3)= T(2,3) + WTMOL * CMASSY * CMASSZ IF(IFLG) 9,10,10 C IF THIS IS A DERIVATIVE RUN, JUST TRANSFORM INERTIAL MATRIX 9 DO 12 I=1,3 12 CMASS(I)=COSXY(1,I)*(T(1,1)*COSXY(1,I)+2.*T(1,2)*COSXY(2,I)+2.* 1T(1,3)*COSXY(3,I))+COSXY(2,I)*(T(2,2)*COSXY(2,I)+2.*T(2,3)*COSXY( 23,I))+T(3,3)*COSXY(3,I)**2-ROT(I) A=DA*CMASS(1) B=DB*CMASS(2) C=DC*CMASS(3) DO 200 N=1,NOATM CMASS(1)=COORD(N,1)-CMASSX CMASS(2)=COORD(N,2)-CMASSY CMASS(3)=COORD(N,3)-CMASSZ DO 200 K=1,3 SUM=0. DO 201 I=1,3 201 SUM=SUM+COSXY(I,K)*CMASS(I) 200 PACOR(N,K)=SUM RETURN C DIAGONALIZE AND ORDER THE 10 N=3 I=0 CALL PRINCE(ROT,T,COSXY) C PRINT OUT ALL RESULTS CMASS(1)=CMASSX CMASS(2)=CMASSY CMASS(3)=CMASSZ DO 60 N=1,NOATM DO 60 I=1,3 60 CMCOR(N,I)=COORD(N,I)-CMASS(I) if(iopt(4).eq.9)itape2=7 if(iopt(4).eq.0)itape2=6 WRITE (ITAPE2,908)CMASSX,CMASSY,CMASSZ WRITE (ITAPE2,904)(NO(J),(CMCOR(J,I),I=1,3),J=1,NOATM) DO 65 N=1,NOATM DO 66 K=1,3 66 CMASS(K)=CMCOR(N,K) DO 65 K =1,3 SUM = 0.0 DO 64 I=1,3 64 SUM=SUM+COSXY(I,K )*CMASS(I) 65 PACOR(N,K)=SUM 70 WRITE (ITAPE2,905)(NO(J),(PACOR(J,I),I=1,3),J=1,NOATM) WRITE (ITAPE2,150)ROT(1),(COSXY(K,1),K=1,3) WRITE (ITAPE2,151)ROT(2),(COSXY(K,2),K=1,3) WRITE (ITAPE2,152)ROT(3),(COSXY(K,3),K=1,3) A = 505376./ROT(1) B = 505376./ROT(2) C = 505376./ROT(3) DA=-A/ROT(1) DB=-B/ROT(2) DC=-C/ROT(3) ASYPAR = (2.0*B - A-C)/(A-C) IF(ASYPAR)76,77,77 76 BPAR=(C-B)/(2.0*A-B-C) GO TO 78 77 BPAR=(A-B)/(2.0*C-B-A) 78 PAR1=B+C PAR2=A+C PAR3=A-C PAR4=B-C PAR5=A-B PAR6=2.0*A-B-C PAR7=2.0*B-A-C PAR8=2.0*C-B-A WRITE (ITAPE2,153)A,PAR1,B,PAR2,C,PAR3,ASYPAR,PAR4,BPAR,PAR5,PAR6, 1PAR7,PAR8 C CALCULATION OF CONSTANTS FOR THE INTERNAL ROTATION PROBLEM IF (ITOP ) 974,11,974 974 DO 40 I=1,NOATM IF(ITOP.EQ.NO(I)) GO TO 41 40 CONTINUE GO TO 11 41 ITOP=I DO 42 I=1,NOATM IF(ITOPA.EQ.NO(I)) GO TO 43 42 CONTINUE GO TO 11 43 ITOPA=I DELCO1 = PACOR(ITOP,1) - PACOR(ITOPA,1) DELCO2 = PACOR(ITOP,2) - PACOR(ITOPA,2) DELCO3 = PACOR(ITOP,3) - PACOR(ITOPA,3) VECLEN=DSQRT(DELCO1**2+DELCO2**2+DELCO3**2) ALAMDA = DELCO1/VECLEN BLAMDA = DELCO2/VECLEN CLAMDA = DELCO3/VECLEN DO 44 I=1,NOATM IF(NA(I).EQ.ITOP.AND. I.NE.ITOPA) GO TO 45 44 CONTINUE RETURN 45 R1=PACOR(I,1)-PACOR(ITOP,1) R2=PACOR(I,2)-PACOR(ITOP,2) R3=PACOR(I,3)-PACOR(ITOP,3) ALFAI=ALAMDA*R1+BLAMDA*R2+CLAMDA*R3 ALFAI=R1**2+R2**2+R3**2-ALFAI**2 ALFAI=3.*WTATOM(MASS(I))*ALFAI ALPHA = ALAMDA*ALFAI/ROT(1) BETA = BLAMDA*ALFAI/ROT(2) GAMMA = CLAMDA*ALFAI/ROT(3) ALFSQ = ALPHA**2 BETSQ = BETA**2 GAMSQ = GAMMA**2 R = 1.0 - ALAMDA*ALPHA-BLAMDA*BETA-CLAMDA*GAMMA FMHZ=505376./(R*ALFAI) FCAL = FMHZ*2.859/29979.0 FALFSQ = FMHZ*ALFSQ FBETSQ = FMHZ*BETSQ FGAMSQ = FMHZ*GAMSQ WRITE(ITAPE2,975)ALFAI,ALAMDA,BLAMDA,CLAMDA,ALPHA,BETA,GAMMA,ALFSQ U ,BETSQ,GAMSQ,R,FMHZ,FCAL,FALFSQ,FBETSQ,FGAMSQ 11 RETURN 150 FORMAT(1H1,25X,68HTHE PRINCIPAL 1ION COEFFICIENTS ARE//28X,7HMOMENTS,28X,27HTRANSFORMATION COEFFICI 2ENTS//50X,1HX,25X,1HY,25X,1HZ//9X,7HI(A) = ,E17.8,9X,E17.8,9X,E17. 38,9X,E17.8) 151 FORMAT(9X,7HI(B) = ,E17.8,9X,E17.8,9X,E17.8,9X,E17.8) 152 FORMAT(9X,7HI(C) = ,E17.8,9X,E17.8,9X,E17.8,9X,E17.8) 153 FORMAT(1H0,23X,28HTHE ROTATIONAL CONSTANTS ARE//35X,4HA = ,F12.3,2 17X,6HB+C = ,F12.3/35X,4HB = ,F12.3,27X,6HA+C = ,F12.3/35X,4HC = ,F 212.3,27X,6HA-C = ,F15.5/31X,8HKAPPA = ,F12.8,27X,6HB-C = ,F15.5/24 3X,15HB(P) OR B(O) = ,F12.8,27X,6HA-B = ,F15.5/75X,9H2A-B-C = ,F15. 45/75X,9H2B-A-C = ,F15.5/75X,9H2C-B-A = ,F15.5) 908 FORMAT( ///6X, 66HTHE CENTER OF MASS COORDINATES IN THE FRA 2ME DEFINED ABOVE ARE X= ,F10.5,5H ,Y= ,F10.5,5H ,Z= ,F10.5////) 904 FORMAT(34X,52HTHE ATOM COORDINATES IN THE CENTER OF MASS FRAME ARE 1//8X,8HATOM NO.,19X,1HX,23X,1HY,23X,1HZ, //(9X,I5,13X,F18 2.9,6X,F18.9,6X,F18.9)) 905 FORMAT(1H0,36X,48HTHE ATOM COORDINATES IN THE PRINCIPLE AXIS FRAME 1//8X,8HATOM NO.,19X,1HA,23X,1HB,23X,1HC, //(9X,I5,13X,F18 2.9,6X,F18.9,6X,F18.9)) 975 FORMAT (1H5,40H CONSTANTS FOR INTERNAL ROTATION PROBLEM/// 5H TOP 1 9HINERTIA =F12.3,/,12H LAMBDA A =, 2 F8.5/12H LAMBDA B = ,F8.5/12H LAMBDA C = ,F8.5/// 9H ALP 3HA = ,F9.6/9H BETA = ,F9.6/9H GAMMA = ,F9.6/// 12H ALPHA**2 = ,E 412.5/12H BETA**2 = ,E12.5/12H GAMMA**2 = ,E12.5/// 55H R = ,F8.5/5H F = ,F13.2, 6 9H MHZ = ,F9.3,15H CALORIES/MOLE/// 15H F*ALPHA**2 = , 7F9.3/15H F*BETA**2 = ,F9.3/15H F*GAMMA**2 = ,F9.3) END SUBROUTINE PRINCE(D,H,U) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) DIMENSION H(3,3),U(3,3),D(3),A(3) DO 100 I=1,3 DO 101 J=1,3 101 U(I,J)=0. D(I)=H(I,I) 100 U(I,I)=1. IFLAG=-2 A(1)=H(2,3) A(2)=H(1,3) A(3)=H(1,2) C FIND LARGEST OFF DIAGONAL ELEMENT PIV=DABS(A(1)) C =DABS(A(2)) S =DABS(A(3)) NPV=1 IPV=2 LPV=3 IF(PIV.GT.C.AND.PIV.GT.S) GO TO 1 IPV=1 IF(S.GT.C) GO TO 103 NPV=2 LPV=3 PIV=C GO TO 1 103 NPV=3 LPV=2 PIV=S 1 I=LPV LPV=IPV IPV=NPV NPV=I DS=D(NPV)-D(LPV) AS=DABS(DS) C CHECK TO SEE IF LARGEST ELEMEMT HAS A SIGNIFICANT ROTATION IF(PIV.LE.1.0E-08*AS) GO TO 7 2 IFLAG=-1 C CHECK FOR SPECIAL CASE OF A VERY BIG PERTURBATION IF(1.0E-10*PIV.LT.AS) GO TO 3 4 S = 0.7071067811865D 00 C=S GO TO 5 C CALCULATE ROTATION 3 T=A(IPV)/DS S=0.25/DSQRT(0.25+T*T) C=DSQRT(0.5+S) S=(T+T)*S/C C ROTATE 5 CC=C*C SS=S*S CS=C*S PIV=(CS+CS)*A(IPV) T=D(LPV) D(LPV)=T*CC+D(NPV)*SS-PIV D(NPV)=T*SS+D(NPV)*CC+PIV A(IPV)=(CC-SS)*A(IPV)-DS*CS T=A(LPV) A(LPV)=S*A(NPV)+C*T A(NPV)=C*A(NPV)-S*T DO 6 K=1,3 T=U(K,LPV) U(K,LPV)=C*T-U(K,NPV)*S 6 U(K,NPV)=S*T+U(K,NPV)*C C FIND THE NEW LARGEST ELEMENT PIV=DABS(A(NPV)) IF(PIV.GT.DABS(A(LPV)) ) GO TO 1 PIV=DABS(A(LPV)) I=LPV LPV=NPV NPV=I GO TO 1 C FIND A ROTATION THAT IS SIGNIFICANT 7 IF(IFLAG.EQ.0) GO TO 10 IFLAG=IFLAG+1 PIV=DABS(A(NPV)) GO TO 1 C ORDER THE EIGENVALUES AND DO THE CORRECT THING TO THE TRANSFORMATI 10 DO 11 I=2,3 II=5-I DO 12 J=2,II IF(D(J).GE.D(J-1) ) GO TO 12 T=D(J) D(J)=D(J-1) D(J-1)=T DO 13 K=1,3 T=U(K,J-1) U(K,J-1)=-U(K,J) 13 U(K,J)=T 12 CONTINUE 11 CONTINUE RETURN END SUBROUTINE DIPOLE(ISIGNL) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) INTEGER*4 ISIGNL LOGICAL NRGTST DIMENSION NPOS(30),NNEG(30),BONDMO(30),DIPOL(145,4),APROJ(30) DIMENSION BPROJ(30),CPROJ(30),BDMIDA(30),BDMIDB(30),BDMIDC(30) COMMON XYZ(50,3),DUM(320),NATM,NO(51),MUD(50),NA(50),NB(50), 1NC(50),NRDN,NZ,IOPT(4) IF (ISIGNL)500,501,107 107 IF (IOPT(1).LT.0) GO TO 101 NOBDDI=IOPT(1) READ (5,102) (NPOS(I),NNEG(I),BONDMO(I),I=1,NOBDDI) 102 FORMAT (2I5,F10.3) DO 103 I=1,NOBDDI ICTPOS=0 ICTNEG=0 DO 106 J=1,NATM IF (NO(J).NE.NPOS(I).OR.ICTPOS.NE.0)GO TO 104 NPOS(I)=J ICTPOS=1 GO TO 105 104 IF (NO(J).NE.NNEG(I).OR.ICTNEG.NE.0)GO TO 105 NNEG(I)=J ICTNEG=1 105 IF (ICTPOS.NE.0.AND.ICTNEG.NE.0)GO TO 103 106 CONTINUE 103 CONTINUE 101 K=1 ENERGY=0. GO TO 110 501 IF (IOPT(1).NE.0)GO TO 502 WRITE (6,503) 503 FORMAT (//52H DIPOLE MOMENT DERIVATIVES CANNOT BE OBTAINED WITHOU, 138HT REQUESTING DIPOLE MOMENTS THEMSELVES) CALL EXIT 502 K=K+1 110 ADIP=0. BDIP=0. CDIP=0. DO 120 I=1,NOBDDI A=XYZ(NPOS(I),1)-XYZ(NNEG(I),1) B=XYZ(NPOS(I),2)-XYZ(NNEG(I),2) C=XYZ(NPOS(I),3)-XYZ(NNEG(I),3) BDLTH=DSQRT(A**2+B**2+C**2) APROJ(I)=BONDMO(I)*A/BDLTH BPROJ(I)=BONDMO(I)*B/BDLTH CPROJ(I)=BONDMO(I)*C/BDLTH ADIP=ADIP+APROJ(I) BDIP=BDIP+BPROJ(I) CDIP=CDIP+CPROJ(I) IFD2=IOPT(2)**2 NRGTST=K.NE.1.OR.IFD2.LE.1 IF (NRGTST) GO TO 120 BDMIDA(I)=XYZ(NPOS(I),1)+XYZ(NNEG(I),1) BDMIDB(I)=XYZ(NPOS(I),2)+XYZ(NNEG(I),2) BDMIDC(I)=XYZ(NPOS(I),3)+XYZ(NNEG(I),3) IF (I.EQ.1) GO TO 120 IMINUS=I-1 DO 130 J=1,IMINUS AONE=(BDMIDA(I)-BDMIDA(J))/2. BONE=(BDMIDB(I)-BDMIDB(J))/2. CONE=(BDMIDC(I)-BDMIDC(J))/2. ONELTH=DSQRT(AONE**2+BONE**2+CONE**2) AONE=AONE/ONELTH BONE=BONE/ONELTH CONE=CONE/ONELTH ENERGY=ENERGY+(APROJ(I)*APROJ(J)+BPROJ(I)*BPROJ(J)+CPROJ(I)* 1CPROJ(J)-3*(APROJ(I)*AONE+BPROJ(I)*BONE+CPROJ(I)*CONE)*(APROJ(J)* 2AONE+BPROJ(J)*BONE+CPROJ(J)*CONE))/(ONELTH*ONELTH*ONELTH) 130 CONTINUE 120 CONTINUE DIPOL(K,2)=ADIP DIPOL(K,3)=BDIP DIPOL(K,4)=CDIP IF (K.EQ.1) GO TO 135 DIPOL(K,2)=1000.*(ADIP-DIPOL(1,2)) DIPOL(K,3)=1000.*(BDIP-DIPOL(1,3)) DIPOL(K,4)=1000.*(CDIP-DIPOL(1,4)) DIPOL(K,1)=(DIPOL(1,2)*DIPOL(K,2)+DIPOL(1,3)*DIPOL(K,3)+ 1 DIPOL(1,4)*DIPOL(K,4))/DIPOL(1,1) RETURN 135 WRITE (6,600)NOBDDI 600 FORMAT (//15X,13HTHE FOLLOWING,I3,28H POLAR BONDS CONTRIBUTE TO T, 116HHE DIPOLE MOMENT//1X9HPOS. ATOM,3X9HNEG. ATOM,5X11HBOND MOMENT, 26X11HA COMPONENT,6X11HB COMPONENT,6X11HC COMPONENT/) DIPOL(1,1)=DSQRT(ADIP**2+BDIP**2+CDIP**2) DO 140 I=1,NOBDDI NP=NO(NPOS(I)) NN=NO(NNEG(I)) WRITE (6,601)NP,NN,BONDMO(I),APROJ(I),BPROJ(I),CPROJ(I) 140 CONTINUE 601 FORMAT (5X,I2,10X,I2,4(10X,F7.3)) WRITE(6,602)DIPOL(1,1),ADIP,BDIP,CDIP 602 FORMAT (//28X34HMOLECULAR DIPOLE MOMENT COMPONENTS//1X, 110HMU(TOTAL)=,F6.3,12X6HMU(A)=,F6.3,12X6HMU(B)=,F6.3,12X6HMU(C)=, 2F6.3) IF (IFD2-4)150,160,170 160 WRITE (6,603) 603 FORMAT (/1X51HTHE TOTAL DIPOLE-DIPOLE INTERACTION ENERGY OF THIS , 160HCONFORMATION HAS BEEN SET EQUAL TO 0.0 KCAL/MOLE AND WILL BE/1X 252HUSED AS A REFERENCE FOR THE FOLLOWING CONFORMATIONS.) REFRNC=ENERGY RETURN 170 ENERGY=14.39*(ENERGY-REFRNC) WRITE (6,604) ENERGY 604 FORMAT (/1X38HTHE DIPOLE-DIPOLE INTERACTION ENERGY =,F6.2,5H KCAL, 147H/MOLE RELATIVE TO THE LAST PRECEDING REFERENCE.) 150 RETURN 500 WRITE (6,605) 605 FORMAT (/1X51HDERIVATIVES OF DIPOLE MOMENT COMPONENTS WITH RESPEC, 122HT TO INPUT COORDINATES/5X34HUNITS ARE MILLI-DEBYE PER DEGREE O 216HR MILLI-ANGSTROM/7X20HDISTANCE DERIVATIVES,22X13HTHETA DERIVAT, 34HIVES,25X15HPHI DERIVATIVES/1X,3HATM,3(3X,3HATM,2X,7HMU(TOT),3X, 45HMU(A),4X,5HMU(B),4X,6HMU(C) )) DO 190 I=1,4 DIPOL(1,I)=0.0 190 CONTINUE KPHI=1 NPHI=0 KDIST=2 KTHETA=1 NTHETA=0 DO 180 I=2,NATM NUMBI=NO(I) NDIST=NO(NA(I)) IF (I-4)191,192,193 191 IF (I.EQ.2)GO TO 195 KDIST=3 KTHETA=4 GO TO 194 192 KDIST=5 KTHETA=6 KPHI=7 GO TO 196 193 KDIST=KDIST+3 KTHETA=KTHETA+3 KPHI=KPHI+3 196 NPHI=NO(NC(I)) 194 NTHETA=NO(NB(I)) 195 WRITE(6,606)NUMBI,NDIST,(DIPOL(KDIST,J),J=1,4),NTHETA, 1 (DIPOL(KTHETA,J),J=1,4),NPHI,(DIPOL(KPHI,J),J=1,4) 606 FORMAT (1X,I2,3(4X,I2,4F9.3)) 180 CONTINUE RETURN END SUBROUTINE ROTOR(A,B,C) IMPLICIT REAL*8 (A-H,O-Z) , INTEGER (I-N) INTEGER TAU DIMENSION XY(50,6) COMMON XZ(50,3),Q(300),WT(20),NATM,NO(51),MASS(204),NTOP,IOPT COMMON/DUMCOM/X1(50,3),X2(50,3),T(3,3), ROT(3),XC(3),XCC(3),XCR(3, C 3),E(3),R(3),RR(3),RINV(3,3),RC(3,3),NX(52) EQUIVALENCE (X1,XY) READ(5,100) NSTEPS,NSIZE,NEG,NPOS,(NX(I),I=1,NTOP) 100 FORMAT(16I5) if(iopt.eq.9)go to 1003 WRITE(6,101) NEG,NPOS,(NX(I),I=1,NTOP) 1003 continue 101 FORMAT(///,50H CALCULATIONS FOR AN INTERNAL ROTOR ALONG AN AXIS F 9HFROM ATOM,I3, 8H TO ATOM ,I3,/29H THE FOLLOWING ATOMS ARE NAME F 20HD AS PART OF THE TOP,/,(20I5)) C FIND AXIS ATOMS W=0. MPOS=1 MNEG=2 DO 10 I=1,NATM W=W+WT(MASS(I)) IF(NPOS.EQ.NO(I)) MPOS=I IF(NEG.EQ.NO(I)) MNEG=I DO 10 J=1,3 X1(I,J)=0. 10 X2(I,J)=0. R(1)=XZ(MPOS,1) R(2)=XZ(MPOS,2) R(3)=XZ(MPOS,3) E(1)=R(1)-XZ(MNEG,1) E(2)=R(2)-XZ(MNEG,2) E(3)=R(3)-XZ(MNEG,3) ER=DSQRT(E(1)**2+E(2)**2+E(3)**2) DO 1 I=1,3 XC(I)=0. E(I)=E(I)/ER DO 1 J=1,3 RINV(I,J)=0. 1 XCR(I,J)=0. C FIND CARTESIAN COORDINATES AS A FUNCTION OF TAU XX=0. DO 2 I=1,NATM DO 3 J=1,NTOP IF(NX(J).EQ.NO(I)) GO TO 4 3 CONTINUE GO TO 2 4 ER=0. DO 5 K=1,3 RR(K)=XZ(I,K)-R(K) 5 ER=ER+E(K)*RR(K) WW=WT(MASS(I)) KK=2 KKK=3 DO 6 K=1,3 X1(I,K)=E(KK)*RR(KKK)-E(KKK)*RR(KK) XC(K)=XC(K)+WW*X1(I,K)/W XX=XX+WW*X1(I,K)**2 X2(I,K)=E(K)*ER-RR(K) KK=KKK 6 KKK=K DO 7 K=1,3 XCR(1,K)=WW*(XZ(I,KK)*X1(I,KKK)-XZ(I,KKK)*X1(I,KK))+XCR(1,K) XCR(2,K)=WW*(XZ(I,KK)*X2(I,KKK)-XZ(I,KKK)*X2(I,KK))+XCR(2,K) KK=KKK 7 KKK=K 2 CONTINUE C FIND TOP INERTIA, X VECTOR DEPENDENCE, ETC. ER=0. ROT(1)=505376./A ROT(2)=505376./B ROT(3)=505376./C KK=2 KKK=3 DO 8 K=1,3 XX=XX-W*XC(K)**2 XCC(K)=E(KK)*XC(KKK)-E(KKK)*XC(KK) R(K)=-XCR(1,K)/ROT(K) ER=ER+R(K)*XCR(1,K) KK=KKK 8 KKK=K XCR(3,1)=E(1)*XX-XCR(1,1) XCR(3,2)=E(2)*XX-XCR(1,2) XCR(3,3)=E(3)*XX-XCR(1,3) GINV=XX+ER TAU=0 RINV(1,1)=A RINV(2,2)=B RINV(3,3)=C C FIND CARTESIAN COORDINATES IN CENTER OF MASS SYSTEM DO 11 I=1,NATM DO 11 J=1,3 X1(I,J)=X1(I,J)-XC(J) 11 X2(I,J)=X2(I,J)-XCC(J) IF(IOPT.EQ.9)GO TO 1005 WRITE(6,102)E,XX,XCR,(NO(I),(XY(I,J),J=1,6),I=1,NATM) 1005 CONTINUE 102 FORMAT(//, 24H AXIS DIRECTION COSINES ,3F10.5,16H TOP INERTIA = 1F10.3,9X,19HX VECTOR COMPONENTS,3(/,80X,3E15.5)//16X,10HCOORDINATE F 23H DEPENDENCE ON SIN(TAU),30X,25H DEPENDENCE ON 1-COS(TAU),/, F 9H ATOM NO.,2(7X,1HA,19X,1HB,19X,1HC,12X),/,(I5,6F20.8)) WRITE(6,104) 104 FORMAT(//5H TAU,3x,9HG INVERSE,4x,1HA,9x,1HB,9x,1HC,8x, F 3HB+C,7x,5hkappa,6x, F 51H ROTATIONAL MATRIX IN FRAME FIXED AXIS C VECTOR) BC=B+C bk=(B+B-A-C)/(A-C) IF(IOPT.GT.0.AND.BK.GT.-0.60)GO TO 12 WRITE(6,103) TAU,GINV,A,B,C,bc,bk,((RINV(I,J),J=I,3),R(I),I=1,3) 12 CONTINUE 103 FORMAT(/,I5,5f10.2,f10.4,8X,3F10.2,E22.5,/,83X,2F10.2,E22.5,/ 93x, f f10.2,e22.5) IF(NSTEPS.LT.1) RETURN C FIND INERTIAL MATRIX, ETC., FOR OTHER VALUES OF TAU DO 20 ISTEP=1,NSTEPS TAU=TAU+NSIZE ANG = 0.0174532925199D 00*TAU CT=1.-DCOS(ANG) ST=DSIN(ANG) DO 30 JSTEP=1,2 DO 21 I=1,3 XC(I)=XCR(1,I)+ST*XCR(2,I)+CT*XCR(3,I) DO 21 J=I,3 21 RC(I,J)=0. DO 22 IN=1,NATM DO 23 I=1,3 23 R(I)=XZ(IN,I)+ST*X1(IN,I)+CT*X2(IN,I) WW=WT(MASS(IN)) DO 22 I=1,3 DO 22 J=I,3 22 RC(I,J)=RC(I,J)-WW*R(I)*R(J) EN=RC(1,1)+RC(2,2)+RC(3,3) RC(1,1)=RC(1,1)-EN RC(2,2)=RC(2,2)-EN RC(3,3)=RC(3,3)-EN RINV(1,1)=RC(2,2)*RC(3,3)-RC(2,3)**2 RINV(2,2)=RC(1,1)*RC(3,3)-RC(1,3)**2 RINV(3,3)=RC(1,1)*RC(2,2)-RC(1,2)**2 RINV(1,2)=RC(1,3)*RC(2,3)-RC(1,2)*RC(3,3) RINV(1,3)=RC(1,2)*RC(2,3)-RC(1,3)*RC(2,2) RINV(2,3)=RC(1,2)*RC(1,3)-RC(2,3)*RC(1,1) EN=RINV(1,1)*RC(1,1)+RINV(1,2)*RC(1,2)+RINV(1,3)*RC(1,3) GINV=XX DO 24 I=1,3 R(I)=0. DO 25 J=1,3 RINV(J,I)=RINV(I,J) 25 R(I)=R(I)-RINV(I,J)*XC(J)/EN 24 GINV=GINV+R(I)*XC(I) CALL PRINCE(ROT,RC,T) EN=505376./EN DO 26 I=1,3 ROT(I)=505376./ROT(I) DO 26 J=I,3 26 RINV(I,J)=RINV(I,J)*EN bd=rot(2)+rot(3) bx=(rot(2)+rot(2)-rot(1)-rot(3))/(rot(1)-rot(3)) IF(IOPT.GT.0.AND.BX.GT.-0.75)GO TO 27 WRITE(6,103) TAU,GINv,rot,bd,bx,((RINV(I,J),J=I,3),R(I),I=1,3) 27 CONTINUE IF(TAU.EQ.180) GO TO 20 TAU=-TAU 30 ST=-ST 20 CONTINUE RETURN END