PROGRAM MOMENT CABCKAP C PARAMETERS A, B, C, AND KAPPA FOR ROTATIONAL SPECTRA C REQUIRES THE FOLLOWING SUBROUTINES--MMPY,MAT,CONST,MPY C COORDINATES ARE GENERATED BY ROTATING AND TRANSPOSING THE C COORDINATE AXES. THE AXIS SYSTEM IS RIGHT HANDED. TRANSPOSING C THE SYSTEM, CORRESPONDING TO PRODUCING THE BOND LENGTH, IS DONE C ALONG THE X AXIS ONLY. A POSITIVE BOND LENGTH TRANSPOSES X TO C X MINUS BOND LENGTH. C ROTATION ABOUT THE AXES IS DONE BY ROTATING THE AXES CLOCKWISE, C FIRST ABOUT THE Z AXIS, SECOND ABOUT THE X AXIS, THIRD ABOUT THE C Y AXIS. CHARACTER*72 TITLE dimension AW(100),CD(3,100),C(3,100),B(3,3),V(3,3),VT(3,3), 1U(3,3),P(3,3),THETA(100),TH(100),PHI(100),PH(100), 2ALPHA(100),AL(100),BONDL(100),ROTT(100,3,3 ),ROX(100,3,3), 3 ROTP(100,3,3), ROTAL(100,3,3), 4 ALX(20), THX(20), PHX(20), BONDX(20), CC1(20), CC2(20), 5 CC3(20), IMCH(20), CMAS(20), RO(100,3,3), CX(3,100), PAR(3) dimension ISTP(100), IST(100), IABCH(20),index(3),ICD(20) OPEN(5,FILE='TAPE5.txt') OPEN(6,FILE='TAPE6.txt') 1 NOCHD=0 ICC1=0 ICC2=0 ICA1=0 ICA2=0 c READ(5,110)TITLE READ(5,100) N, IFAB, NOCH, DELTA WRITE(6,100) N, IFAB, NOCH, DELTA C IFAB=1 IF ANGLES AND LENGTHS ARE FED IN, OTHERWISE 0 C NOCH= NUMBER OF RECALCULATIONS TO BE DONE IF (N.LE.0) GO TO 50 write(6,*) 'IFAB=',ifab IF (IFAB.EQ.0) GO TO 2 C DEPARTURE POINT WHEN MOLECULE IS CONSTRUCTED DO 210 I=1,N 210 READ(5,120) ISTP(I), AL(I), TH(I), PH(I), BONDL(I), AW(I) 211 DO 215 I=1,N ALPHA(I)=AL(I)*0.017453293 THETA(I)=TH(I)*0.017453293 215 PHI(I)=PH(I)*0.017453293 CALL CONST(ALPHA,1,2,N,ROTAL) CALL CONST (THETA,2,3,N,ROTT) CALL CONST (PHI,3,1,N,ROTP) DO 220 I=1,N DO 220 J=1,3 DO 220 L=1,3 ROX(I,J,L)=0. DO 220 K=1,3 220 ROX(I,J,L)=ROX(I,J,L)+ROTP(I,J,K)*ROTT(I,K,L) DO 225 I=1,N DO 225 J=1,3 DO 225 L=1,3 RO(I,J,L)=0. DO 225 K=1,3 225 RO(I,J,L)=RO(I,J,L)+ROX(I,J,K)*ROTAL(I,K,L) IA=0 IB=0 230 IF (IA.GE.N) GO TO 4 IA=IA+1 IA1=IA-1 IF(IA1.LE.0) GO TO 280 IF (ISTP(IA)-IB) 245,245,270 C BEGIN TO RETRACE WHEN PREVIOUS STARTING POINT CALLED FOR 245 ISTX=IST(IB) DO 260 I=1,IA1 260 CD(1,I)=CD(1,I)+BONDL(ISTX) CALL MPY (RO,CD,1,ISTX,IA1,CD) IB=IB-1 IF (ISTP(IA)-IB) 245,245,270 C END OF RETRACE TO PREVIOUS STARTING POINT C NOW ADDITION OF IA-TH ATOM 270 CALL MPY(RO,CD,0,IA,IA1,CD) DO 275 I=1,IA1 275 CD(1,I)=CD(1,I)- BONDL(IA) 280 DO 285 I=1,3 285 CD(I,IA)=0. IB=IB+1 IST(IB)=IA GO TO 230 C END OF CALCULATION OF COORDINATES C START HERE WHEN COORDINATES ARE FED IN 2 DO 3 I=1,N 3 READ(5,102) AW(I),CD(1,I), CD(2,I),CD(3,I) WRITE(6,102) AW(I),CD(1,I), CD(2,I),CD(3,I) 4 WT=0. c WRITE(6,110)TITLE DO 5 I=1,N 5 WT=WT+AW(I) DO 10 J=1,3 CM=0. DO 8 I=1,N 8 CM=CM+CD(J,I)*AW(I) DO 9 I=1,N 9 C(J,I)=CD(J,I)-CM/WT 10 CONTINUE CALL MAT(0,B) DO 15 I=1,N B(1,1)=B(1,1)+AW(I)*C(2,I)**2+AW(I)*C(3,I)**2 B(2,2)=B(2,2)+AW(I)*C(1,I)**2+AW(I)*C(3,I)**2 B(3,3)=B(3,3)+AW(I)*C(1,I)**2+AW(I)*C(2,I)**2 DO 15 K=1,3 DO 15 L=1,3 IF (K.EQ.L) GO TO 15 B(K,L)=B(K,L)-AW(I)*C(K,I)*C(L,I) 15 CONTINUE C DIAGONALIZE INERTIA TENSOR S=(AMAX1(B(1,1),B(2,2),B(3,3)))*1.E-8 CALL MAT(1,U) 19 OFF=0. DO 29 I=1,2 II=I+1 DO 29 J=II,3 IF ( ABS(B(I,J)).LE.OFF) GO TO 29 OFF=ABS(B(I,J)) IO=I JO=J 29 CONTINUE IF (OFF.LE.S) GO TO 32 CALL MAT(1,V) CALL MAT(1,VT) DD=(ABS (B(IO,IO)-B(JO,JO)))*10. IF (DD.GE.B(IO,JO)) GO TO 31 ANG = 0.78539816 GO TO 1132 31 ARG =(2.*B(IO,JO))/(B(IO,IO)-B(JO,JO)) ANG=(ATAN(ARG))/2. 1132 SANG = SIN(ANG) CANG =COS(ANG) V(IO,JO)=-SANG V(JO,IO)=SANG V(IO,IO)=CANG V(JO,JO)=CANG VT(IO,JO)=SANG VT(JO,IO)=-SANG VT(IO,IO)=CANG VT(JO,JO)=CANG CALL MMPY(B,V,P) CALL MMPY (VT,P,B) CALL MMPY(VT,U,U) GO TO 19 32 DO 33 I=1,3 DO 33 J=1,N CX(I,J)=0. DO 33 K=1,3 33 CX(I,J)=CX(I,J)+U(I,K)*C(K,J) DO 34 I=1,3 DO 34 J=1,N 34 C(I,J) =CX(I,J) DO 1334 K=1,3 BIGER =0. DO 1234 I=1,3 IF(BIGER.GE.B(I,I)) GO TO 1234 BIGER =B(I,I) ISAV=I INDEX(K)=I PAR(K)=BIGER 1234 CONTINUE 1334 B(ISAV,ISAV)=0. AI=PAR(3) BI=PAR(2) CI=PAR(1) AAM=BI+CI-AI BBM=AI+CI-BI CCM=AI+BI-CI PA=505379.1/AI PB=505379.1/BI PC=505379.1/(CI+DELTA) MM1=INDEX(3) MM2=INDEX(2) MM3=INDEX(1) PD=(PA-PC)/2.0 PE=(PA+PC)/2.0 ASYM=(2.*PB-PA-PC)/(PA-PC) BANDC = PB+PC IF (IFAB.EQ.0) GO TO 400 C WRITE(6,112) WRITE(6,121) C WRITE(6,111) WRITE(6,122) C WRITE(6,111) DO 501 I=1,N 501 WRITE(6,125) ISTP(I), AL(I), TH(I), PH(I), BONDL(I), AW(I) 400 WRITE(6,112) WRITE(6,109) WT C WRITE(6,111) WRITE(6,105) N C WRITE(6,111) WRITE(6,4000) 4000 FORMAT (21X,1HA,9X,1HB,9X,1HC,14X,1HX,9X,1HY,9X,1HZ) C DO 45 I=1,N DO 45 I=1,N C 1CD(3,I) 45 WRITE(6,103) AW(I),C(MM1,I),C(MM2,I),C(MM3,I),CD(1,I),CD(2,I), 1CD(3,I) c for the ester CH3 top axis is about the C1-C5 direction (low barrier) c for the ester CH3 the other top axis is about the O13-C14 direction (high barrier) c for the carbamate CH3 top axis is about the C7-O1 direction c for the eaac1 CH3 top axis is about the C3-C4 direction c for the eaacclean CH3 top axis is about the C1-C4 direction c for the Phe-Ala-Met-ester low barrier C1-H1-H2-H3 c for the Phe-Ala-Met-ester high barrier C10-H11-H12-H13 c A2=C(MM1,13) A2=C(MM1,5) c A1=C(MM1,14) A1=C(MM1,1) c A1=C(MM1,10) c B2=C(MM2,13) B2=C(MM2,5) c B1=C(MM2,14) B1=C(MM2,1) c B1=C(MM2,10) c C2=C(MM3,13) C2=C(MM3,5) c C1=C(MM3,14) C1=C(MM3,1) c C1=C(MM3,10) c for the ester we now calculate the center of mass of CH3 cgx=(1/3.)*(C(MM1,2)+C(MM1,3)+C(MM1,4)) c cgx=(1/3.)*(C(MM1,11)+C(MM1,12)+C(MM1,13)) c cgx=(1/3.)*(C(MM1,15)+C(MM1,16)+C(MM1,17)) A2=cgx cgy=(1/3.)*(C(MM2,2)+C(MM2,3)+C(MM2,4)) c cgy=(1/3.)*(C(MM2,11)+C(MM2,12)+C(MM2,13)) c cgy=(1/3.)*(C(MM2,15)+C(MM2,16)+C(MM2,17)) B2=cgy cgz=(1/3.)*(C(MM3,2)+C(MM3,3)+C(MM3,4)) c cgz=(1/3.)*(C(MM3,11)+C(MM3,12)+C(MM3,13)) c cgz=(1/3.)*(C(MM3,15)+C(MM3,16)+C(MM3,17)) C2=cgz c for the acetic acid CH3 top axis is about the C2-C5 direction VR=(A2-A1)**2+(B2-B1)**2+(C2-C1)**2 Vr=sqrt(vr) print*,'VR=',VR,'A2=',A2,'A1=',A1,'B2=',B2,'B1=',B1 1,'C2=',C2,'C1=',C1 rlambA=(A2-A1)/VR print*,'rlamba=',rlamba rlambA=(rlambA*3.2)/AI rlambB=(B2-B1)/VR print*,'rlambb=',rlambb rlambB=(rlambB*3.2)/BI rlambC=(C2-C1)/VR print*,'rlambc=',rlambc rlambC=(rlambC*3.2)/CI print*,'rhoA=',rlambA,'rhoB=',rlambB,'rhoC=',rlambC rho=rlambA**2+rlambB**2+rlambC**2 rho=sqrt(rho) print*,'rho=',rho rr=1.-((rlambA**2)*3.2)/AI - ((rlambB**2)*3.2)/BI c attention the sign of the angle for the two CH3 tops in ester is c opposite: change the sign of rlambda! DAB=(-rlambA*rlambB*(PA-PB))/(rho*rho) DAC=(-rlambA*rlambC*(PA-PC))/(rho*rho) DBC=(rlambB*rlambC*(PB-PC))/(rho*rho) c PAB=505379.1/DAB PF=(1./rr)*505379.1 PF=PF/3.2 PFCM=PF/29979.2458 DABM=DAB/29979.2458 DACM=DAC/29979.2458 DBCM=DBC/29979.2458 c dada=DAB/(2.*(PA-PB)) c print*,'dada=DAB/(2.*(PA-PB))',dada c beta1=sin(dada*0.017453293) c beta2=cos(dada*0.017453293) c beta2=tan(dada*0.017453293) c beta=(atan(dada))*0.5 c print*,'beta en radians=',beta c betad=(beta)/0.01745293 c print*,'beta en degree=',betad c print*,'beta=',beta c ARAM=cos(beta)*PA+sin(beta)*PB c BRAM=-sin(beta)*PA+cos(beta)*PB print*,'rr=',rr,'DAB=',DAB,'PF=',PF,'PFCM=',PFCM print*,'DAB in cm-1=',DABM c print*,'betad=',betad,'ARAM=',ARAM,'BRAM=',BRAM, print*,'PA=',PA 1,'PB=',PB,'PC=',PC CONV=29979.2458 c ARAM2=ARAM/CONV c BRAM2=BRAM/CONV PA2=PA/CONV PB2=PB/CONV PC2=PC/CONV c print*,'ARAM2=',ARAM2,'BRAM2=',BRAM2 print*,'PAM A cm-1=',PA2 print*,'PAM B cm-1=',PB2 print*,'PAM C cm-1=',PC2 WRITE(6,4001) 4001 FORMAT(30H0PRINCIPAL AXIS TRANSFORMATION//10X,1HX,15X,1HY,15X,1HZ) WRITE(6,4002)(U(MM1,I),I=1,3),(U(MM2,I),I=1,3),(U(MM3,I),I=1,3) 4002 FORMAT (' A',3E16.8/' B',3E16.8/' C',3E16.8) WRITE(6,4003) AI,BI,CI,AAM,BBM,CCM,DELTA 4003 FORMAT(' MOMENTS OF INERTIA (AMU*A**2)'/'0I(A) =',F12.7,5X,'I(B) 1=',F12.7,5X,'I(C) =',F12.7/'0I(B) + I(C) - I(A) =',F12.7/' I(A) + 2I(C)- I(B)=',F12.7/' I(A)+I(B)-I(C) =',F12.7,'DELTA = ',F10.7) WRITE(6,4004) 4004 FORMAT(59H0ROTATIONAL CONSTANTS IN MC/S (CONVERSION FACTOR = 5053 176)) WRITE(6,111) WRITE(6,106) PA, PB, PC WRITE(6,111) WRITE(6,107) PD, PE, ASYM, BANDC IF (ICA2-ICA1) 48,48,350 48 IF(ICC2-ICC1) 300,300,387 300 IF (NOCH-NOCHD) 1,1,310 C ENTER CHANGE OF MASS, ANGLE, LENGTH, OR COORDINATE HERE 310 NOCHD=NOCHD+1 ICC1=0 ICC2=0 ICA1=0 ICA2=0 READ(5,130) NCH, MSOC C MSOC = 0 FOR MASS CHANGE, 1 FOR COORDINATE, ANGLE, OR BOND C ENTER MASS CHANGE HERE C IMCH IS THE NUMBER OF THE ATOM FOR WHICH THE MASS IS TO BE CHANGED IF (MSOC.NE.0) GO TO 330 DO 320 I=1,NCH 320 READ(5,132) IMCH(I),CMAS(I) DO 325 I=1,NCH IMCH1=IMCH(I) 325 AW(IMCH1)=CMAS(I) GO TO 4 330 IF (IFAB.EQ.0) GO TO 380 C ENTER BOND LENGTH OR ANGLE CHANGE HERE C IABCH IS NUMBER OF ATOM FOR LENGTH OR ANGLE CHANGE C ICA2 IS THE NUMBER OF CYCLES TO MAKE WITH THIS CHANGE DO 340 I=1,NCH 340 READ(5,134) IABCH(I), ALX(I), THX(I), PHX(I), BONDX(I), ICA2 350 ICA1=ICA1+1 DO 360 I=1,NCH IABC=IABCH(I) AL(IABC)=AL(IABC)+ALX(I) TH(IABC)=TH(IABC)+THX(I) PH(IABC)=PH(IABC)+PHX(I) 360 BONDL(IABC)=BONDL(IABC)+BONDX(I) GO TO 211 C ENTER COORDINATE CHANGE HERE C ICD IS NUMBER OF ATOM WHOSE COORDINATE IS CHANGED 380 DO 385 I=1,NCH 385 READ(5,136) ICD(I),CC1(I),CC2(I),CC3(I),ICC2 387 ICC1=ICC1+1 DO 390 I=1,NCH ICD1=ICD(I) CD(1,ICD1)=CD(1,ICD1)+CC1(I) CD(2,ICD1)=CD(2,ICD1)+CC2(I) 390 CD(3,ICD1)=CD(3,ICD1)+CC3(I) GO TO 4 50 CONTINUE 100 FORMAT (I2,2X,I1,2X,I3,F10.7) 101 FORMAT (I2) 102 FORMAT (F10.7, 3F10.7) 103 FORMAT (F10.5,5X,3F10.6,5X,3F10.6) 105 FORMAT (I3,'ATOMS TRANSPOSED COORDINATES 1ORIGINAL COORDINATES ') 106 FORMAT(5X,'A = ',F11.4,2X,'B = ',F10.4,2X,'C = ',F10.4) 107 FORMAT(5X,'(A-C)/2= ',F11.4,2X,'(A+C)/2 = ',F11.4,2X,'KAPPA = ', 1F11.8,'(B+C) = ',F11.4) 109 FORMAT (1X,19HMOLECULAR WEIGHT = ,F10.5) 110 FORMAT (A72) 111 FORMAT (/) 112 FORMAT (//) 120 FORMAT (I2,8X,3F10.5,F10.7,F10.7) 121 FORMAT(100H0------------------------------------------------------ 1---------------------------------------------) 122 FORMAT (19X,5HALPHA,10X,5HTHETA,11X,3HPHI,9X,11HBOND LENGTH, 15X,7HAT. WT.) 125 FORMAT (I3,8X,3F15.5,F15.7,F15.7) 130 FORMAT (I2,2X,I1) 132 FORMAT (I2,8X,F10.5) 134 FORMAT (I2,8X,3F10.5,F10.7,2X,I3) 136 FORMAT (I2,8X,3F10.7,2X,I3) STOP END SUBROUTINE MMPY(AMP,BMP,CMP) DIMENSION AMP(3,3),BMP(3,3),CMP(3,3),DMP(3,3) DO 10 I=1,3 DO 10 J=1,3 DMP(I,J)=0. DO 10 K=1,3 10 DMP(I,J)=DMP(I,J)+AMP(I,K)*BMP(K,J) DO 20 I=1,3 DO 20 J=1,3 20 CMP(I,J)=DMP(I,J) RETURN END SUBROUTINE MAT(M,BAT) DIMENSION BAT(3,3) DO 10 I=1,3 DO 10 J=1,3 10 BAT(I,J)=0. IF (M.EQ.0) GO TO 20 DO 19 I=1,3 19 BAT(I,I)=1. 20 RETURN END SUBROUTINE CONST (ANGLE,NA,NB,N,ROTAT) DIMENSION ROTAT(100,3,3), ANGLE (100) DO 30 IY=1,N DO 10 I=1,3 DO 10 J=1,3 10 ROTAT(IY,I,J)=0. DO 20 I=1,3 20 ROTAT (IY,I,I)=1. ARGU=ANGLE(IY) ROTAT(IY,NA,NA)=COS (ARGU) ROTAT(IY,NB,NB)=ROTAT(IY,NA,NA) ROTAT(IY,NA,NB)=-SIN(ARGU) ROTAT(IY,NB,NA)=-ROTAT(IY,NA,NB) 30 CONTINUE RETURN END SUBROUTINE MPY (ROTA,CXD,INVERS,ISTY,IX,CYD) DIMENSION ROTA(100,3,3), ROTAX(3,100), CXD(3,100), CYD(3,100) IF (INVERS)10,50,10 10 DO 20 I=1,3 DO 20 J=1,IX ROTAX(I,J)=0. DO 20 K=1,3 20 ROTAX(I,J)=ROTAX(I,J)+ROTA(ISTY,K,I)*CXD(K,J) GO TO 100 50 DO 60 I=1,3 DO 60 J=1,IX ROTAX(I,J)=0. DO 60 K=1,3 60 ROTAX(I,J)=ROTAX(I,J)+ROTA(ISTY,I,K)*CXD(K,J) 100 DO 110 I=1,3 DO 110 J=1,IX 110 CYD(I,J)=ROTAX(I,J) RETURN END