C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C VIBCA - VIBRATIONAL CALCULATIONS PROGRAM C C Developed from P.W.FOWLER's version of VIBCA of 6/1984, with additions C to calculate various matrices of the vibrational problem, PED, C vibrational amplitudes, Coriolis coefficients, quartic cd constants, C inertial defects. C C Companion programs are: C C - VECTOR for displaying normal coordinate displacement vectors C - FCONV for force constant conversions: internal <-> symmetry and C from GAMESS output to VIBCA input C C ver. 11.5.2018 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c 14.06.95: allowance for global scaling factor of force constants c 22.06.95: PARAMETER dimensioning c 10.04.96: disk output of S-reduction quartic cd. constants c 23.04.96: calculation of tau.aaac c 20.06.96: output of force scaling factor and more about units c 23.06.96: elements of H.J and h.3 sextic constants c 19.07.96: correction of external rotation matrix input c 27.03.99: debugging of contributions to delJ + of signs in rotation matrix C 20.10.06: some modifications to output C 27.12.06: MAXAT=100 + related changes in output FORMATs C 8.01.07: Upgrade from the simplified PED calculation which gave erratic C results for lowest frequency modes C 20.05.07: Calculation of vibrational amplitudes C 11.05.18: gfortran compatibility C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C This program takes definitions of internal coordinates and their C associated force constants as data and calculates normal mode frequencies, C eigenvectors (normal coordinate displacement vectors), centrifugal C distortion constants, Coriolis coupling coefficients and harmonic C contributions to moments of inertia. Matrices B, G.int, L.int, L.sym C can also be printed (U matrix is necessary for L.sym). C C Calculation of the normal modes proceeds from internal coordinates via C mass-weighted cartesian coordinates avoiding the definition of symmetry C coordinates and the problem of redundant coordinates. C C see: C C [GWINN] W.D.Gwinn - J.Chem.Phys. 55,477(1971) C C and for calculation of various properties: C C [ALWATS] M.R.Aliev and J.K.G.Watson - J.Mol.Spectrosc. 61,29-52(1976) C [CALIFANO] S.Califano, "Vibrational States", Wiley, London (1976) C [CYVIN] S.J.Cyvin, "Molecular Vibrations and Mean Square Amplitudes", C Elsevier, Amsterdam (1968) C [STOLEVIK] R.Stolevik,H.M.Seip,S.J.Cyvin, Chem.Phys.Lett. 15,263-265(1972) C [WATSCD] J.K.G.Watson - J.Mol.Spectrosc. 65,123-133(1977) C [WATSON] J.K.G.Watson in "Vibrational Spectra and Structure" (J.Durig, C Ed.), Vol.6, Elsevier, Amsterdam. C_____________________________________________________________________________ C C Structure of the input data set: FORMAT C C=============================================== C LINE.1: HDG (20A4) C LINE.2: NAT NFK NLBE IFVCT IPRFUL FSCALE (free format) C LINE.3: AT JJ KK R THETA PHI LL NAX ALFA (1X,A4,2I5,3F10.4,2I5,F10.4) C . C . (atomic coordinate cards - NAT+2*NLBE of these) C . C LINE.I to I+2: optional external 3x3 cart.->princ. rotation matrix C LINE.K: FK(1), FK(2),...,FK(10) (free format - 10 per line) C . . C . .........FK(NFK) C LINE.L: ICC1,LIST(1,1),...,LIST(1,4),INDEX,ICC2,LIST(1,2),...,LIST(4,2) C . . C . . (11I5) C LINE.Ln: 0 C C These lines can optionally be followed by definition of the U matrix C and/or externally supplied moments of inertia and/or experimental C centrifugal distortion constants (in this order) C C LINE U1: U MATRIX (A, capitals only) C LINE U2: row of U matrix containing NINT (free format) C . values, altogether NMOD such rows C . C C LINE I1: MOMENTS OF INERTIA (A, capitals only) C LINE I+1: I.a I.b I.c (free format) C C LINE C1: CENTRIFUGAL DISTORTION CONSTANTS (A, capitals only) C LINE C+1: IASYM D.J D.JK D.K d.J d.K (free format) C C LAST LINE: END OF DATA (A, capitals only) C=============================================== C C The last line terminates the data set, in fact any alphanumeric that C is not one of the header lines is acceptable. C C HDG - TITLE C NAT - NUMBER OF ATOMS IN THE MOLECULE, THIS CAN BE SET TO A NEGATIVE C VALUE IN WHICH CASE EXTERNALLY SUPPLIED ROTATION MATRIX CAN BE C USED TO PRESERVE THE HANDEDNESS OF CARTESIAN SYSTEM ON TRANSITION C TO PRINCIPAL AXES - this is crucial for both reproduction of C literature calculations and of results from GAMESS. C THE 3x3 MATRIX ROTATION MATRIX SHOULD BE PLACED BETWEEN C ATOMIC DEFINITIONS AND FORCE CONSTANT VALUES C NFK - NUMBER OF FORCE CONSTANT VALUES TO BE SPECIFIED C NLBE - NUMBER OF LINEAR BENDS IN THE MOLECULE C (COLLINEAR THREE-ATOM SUBSECTIONS) C IFVCT - EIGENVECTOR SELECTOR, SET TO 1 IF EIGENVECTORS (NORMAL COORDINATES C IN CARTESIANS) ARE TO BE CALCULATED AND OUTPUT, SET TO 0 IF C EIGENVECTORS ARE NOT TO BE PRINTED C IPRFUL - B,G,L, U,L.sym (if U supplied) and Zeta MATRICES TO BE C PRINTED IN BRIEF (=1), IN FULL (=2) OR NOT AT ALL (=0) C FSCALE - force field scaling factor for use with ab-initio force fields, C normally to be set to 1.0 C IASYM - TO BE SET TO -1 FOR PROLATE TYPE CONSTANTS OR TO 1 FOR OBLATE C CONSTANTS C C ATOMIC COORDINATES CAN BE EITHER LOCAL POLARS OR CARTESIANS, C BOTH BEING IN A RIGHT-HANDED AXIS SYSTEM, EACH LINE CONTAINING THE C FOLLOWING: C C AT - SPECIFIES AN ISOTOPE BY ITS MASS NUMBER AND CHEMICAL SYMBOL: if C this isotope is not recognised (see declarations in GEOMAS) then C follow this line with a line containing the mass of this atom C JJ - NUMBER OF ATOM SPECIFYING ORIGIN FOR LOCAL POLAR C COORDINATES OF ATOM I . (JJ=0 FOR ATOM 1 AND FOR CARTESIANS) C KK - NUMBER OF POLAR AXIS (X=1,Y=2,Z=3) ,KK=0 FOR CARTESIANS. C R,THETA,PHI - LOCAL POLAR COORDINATES WITH ORIGIN AT ATOM JJ, C OR CARTESIAN COORDINATES X,Y AND Z. C LL,NAX,ALFA - DESCRIBE A ROTATION OF AXES (SEE SUBROUTINE ROTCS) C C FK - FORCE CONSTANT VALUES (NFK IN NUMBER), TO BE READ WITH TEN VALUES C PLACED IN EACH LINE in units of mdyn/A, mdyn/rad, mdyn A/rad^2 for C stretch, stretch-bend, bend resp. C C LINE.L ONWARDS CONTAIN DEFINITIONS OF INTERNAL COORDINATES AND OF THE C FORCE FIELD, EACH LINE CONTAINS A DEFINITION OF COORDINATES CONTRIBUTING C A TERM TO THE POTENTIAL ENERGY: C C ICC1,ICC2 - CODES SPECIFYING TYPES OF INTERNAL COORDS. 1 AND 2: C 1 = STRETCHING C 2 = BENDING C 3 = OUT-OF-PLANE BENDING C 4 = LINEAR BENDING C 5 = GENERALIZED TORSION ABOUT A BOND C 6 = SIMPLE FOUR ATOM TORSION C LIST - AN ARRAY CONTAINING ATOM NUMBERS DEFINING EACH C INTERNAL COORDINATE C INDEX - SPECIFIES FORCE CONSTANT VALUE FROM ARRAY FK TO BE C ASSOCIATED WITH THE COORDINATES IN THE ENERGY CONTRIBUTION C C IT IS IMPORTANT THAT ALL DIAGONAL TERMS BE SPECIFIED FIRST (otherwise C identification of internal coordinates and the PED calculation may fail) C C I.a I.b I.c - these are optional external values of moments of inertia C (u A**2) to supplant those calculated from the supplied C structures. Internal moments of inertia are used if C the 'MOMENTS OF INERTIA' header line is not encountered C C IASYM D.J D.JK D.K d.J d.K - these are optional experimental values of C quartic centrifugal distortion constants for assessing the C agreement with the calculated values C C----------------------------------------------------------------------------- C C O M P I L A T I O N: C----------------------------------------------------------------------------- C C If you want to compile this program then please note that it may be C dependent on the use of 'static memory allocation'. Please read the following C notes copied from the PROSPE webpage (and also use the data files from C the webpage for testing): C C Most contemporary compilers no longer use default static memory allocation C that preserves values of variables previously set in some subroutine. Many C programs assume this to be the case so that you need to use an appropriate C keyword to enforce static allocation. For example on the f77 compiler on C SGI machines this is the -static option, whereas with Intel Visual Fortran C you have to use option -Qsave. Note that with some compilers optimization C options used not to be safe. With good compilers this is no longer the C case, but if problems crop up it is still good to check by disabling C optimization. C C Command line example for gfortran: C gfortran -fno-automatic vibca.for -o vibca C C Command line example for Intel Visual Fortran 11: C ifort -nopdbfile -nodebug -traceback -O3 -Qsave -ccdefault:fortran C -fpscomp:filesfromcmd vibca.for C C____________________________________________________________________________ C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) DIMENSION HDG(20) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) CHARACTER*30 FILNAM,FILNBM C write(*,'(1x//)') WRITE(*,300) 300 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | V I B C A - VIBrational CAlculations ', * ' using mass-weighted',T79,'|'/ * ' | Cartesian coordinates',T79,'|'/ * ' |',76(1H_),'|'/' version 11.V.2018',T64,'Zbigniew KISIEL'/ * T63,'(Patrick FOWLER)'/) 301 WRITE(*,302)' Name of data file: ' 302 FORMAT(1X/1X,A,$) READ(*,'(A)',ERR=301)FILNAM OPEN(3,FILE=FILNAM,ERR=301,status='old') 303 WRITE(*,302)'Name of output file: ' READ(*,'(A)',ERR=303)FILNBM OPEN(4,FILE=FILNBM,STATUS='UNKNOWN',ERR=303) IR=3 IW=4 write(iw,300) C C C...Read first two lines of data set C READ(IR,101,ERR=301) (HDG(K),K=1,20) WRITE(IW,201) (HDG(K),K=1,20) READ(IR,*) NAT,NFK,NLBE,IFVCT,IPRFUL,FSCALE INPROT=0 IF(NAT.LT.0)THEN NAT=-NAT INPROT=1 ENDIF C C...Read atomic coordinates and convert to principal coordinates C CALL GEOMAS(NAT,NLBE,INPROT) ----- C C...Set up and solve the vibrational problem C CALL VIBFR(NAT,NFK,IFVCT,NMOD,IPRFUL,FSCALE) ----- C C...Calculate cd. constants, zetas and harmonic contributions to moments of C inertia C CALL VRCON(NAT,NMOD,IPRFUL) ----- C WRITE(IW,200) 101 FORMAT(20A4) 102 FORMAT(6I5) 201 FORMAT(1x,77('=')/1x,20A4/1x,77('=')/) 200 FORMAT(1X,77('_')/) C STOP END C C_____________________________________________________________________________ C SUBROUTINE GEOMAS(NAT,NLBE,INPROT) C C GEOMAS HANDLES MOLECULAR GEOMETRY AND ATOMIC MASSES reading data lines of C the form: C C AT JJ KK R THETA PHI LL NAX ALFA (1X,A4,2I5,3F10.4,2I5,F10.4) C C If the isotope definition AT is not included in the recognition table then C this line should be followed by an extra line containing the atomic mass C of that atom. Currently handled isotopes are listed in the DATA EL// C statement. The '*' designation is for dummy atoms involved in linear C angle bends. C C NAT - the number of atoms C NLBE - the number of linear angle bend pairs C INPROT - if set to 1 then external coordinate rotation matrix C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) c CHARACTER POLE(4)*1,EL(28)*4,AT*4 c DIMENSION S(3,3),XP(3),XT(3),CM(3),Z(29) CHARACTER LABEL(3) DATA LABEL /'a','b','c'/ DATA ZERO/0.0D0/,ONE/1.0D0/,RAD/57.295779513082326D0/ DATA EL/' 1H',' 2H',' 3H',' 12C',' 13C',' 14C',' 14N',' 15N', * ' 16O',' 17O',' 18O',' 19F',' 32S',' 34S','35CL','37CL','79BR', * '81BR','127I',' 6LI',' 7LI',' 9BE',' 10B',' 11B','28SI',' 31P', * ' *',' '/ DATA Z/1.00782519D0,2.0141022D0,3.016050D0,12.0D0,13.0033544D0, * 14.003242D0,14.003074D0,15.000107D0,15.994915D0,16.99913D0, * 17.99916D0,18.998405D0,31.97207D0,33.96786D0,34.96885D0, * 36.9659D0,78.9184D0,80.9163D0,126.9043D0,6.01513D0,7.01601D0, * 9.01219D0,10.01294D0,11.009305D0,27.97693D0,30.973765D0, * 3*0.0D0/ DATA POLE/' ','X','Y','Z'/ C C...NATPA - THE NUMBER OF REAL ATOMS PLUS DUMMY ATOMS C NATPA=NAT+NLBE+NLBE DO 1 K=1,NATPA DO 1 KK=1,3 X(KK,K)=ZERO 1 CONTINUE DO 2 I=1,3 DO 3 J=1,3 S(I,J)=ZERO 3 CONTINUE S(I,I)=ONE 2 CONTINUE C C...READ ATOMIC COORDINATES C WRITE(IW,201) DO 4 II=1,NATPA 5 READ(IR,101) AT,JJ,KK,R,THETA,PHI,LL,NAX,ALFA I=II IZ=29 DO 22 M=1,28 IF(EL(M).EQ.AT) IZ=M 22 CONTINUE IF(IZ.EQ.28) I=0 WRITE(IW,202)I,AT,JJ,POLE(KK+1),R,THETA,PHI,LL,NAX,ALFA IF(LL)6,6,7 7 CALL ROTCS(LL,NAX,ALFA,S) ----- IF(I)5,5,6 6 IF(KK)8,8,9 8 X(1,I)=R X(2,I)=THETA X(3,I)=PHI GO TO 23 9 K=MOD(KK,3)+1 L=MOD(KK+1,3)+1 THETA=THETA/RAD PHI=PHI/RAD XP(K)=R*DSIN(THETA) *DCOS(PHI) XP(L)=R*DSIN(THETA)*DSIN(PHI) XP(KK)=R*DCOS(THETA) C C...XP - CARTESIAN COORDINATES IN UNROTATED AXIS SYSTEM WITH C ORIGIN AT ATOM JJ C DO 11 J=1,3 XT(J)=ZERO DO 11 K=1,3 XT(J)=XT(J)+S(J,K)*XP(K) 11 CONTINUE C C...S - MATRIX SPECIFYING AN ORTHOGONAL TRANSFORMATION C XT - CARTESIAN COORDINATES IN ROTATED AXIS SYSTEM C WITH ORIGIN AT ATOM JJ C IF(JJ)12,12,14 12 DO 13 J=1,3 X(J,I)=XT(J) 13 CONTINUE GO TO 23 14 DO 15 J=1,3 X(J,I)=XT(J)+X(J,JJ) 15 CONTINUE 23 IF(IZ.EQ.29) READ(IR,*) Z(IZ) WT(II)=Z(IZ) 4 CONTINUE C C...X - CARTESIANS WITH ORIGIN AT ATOM I C WT - THE ARRAY OF ATOMIC MASSES (AMU) C WM=0 DO 17 I=1,NAT WM=WM+WT(I) 17 CONTINUE C C...WM - THE MOLECULAR WEIGHT C WRITE(IW,203) DO 18 K=1,NATPA WRITE(IW,204) K,(X(KK,K),KK=1,3),WT(K) 18 CONTINUE C C...CALCULATE COORDINATES OF CENTRE OF MASS C DO 19 M=1,3 CM(M) = ZERO DO 20 K=1,NAT CM(M)=CM(M) + WT(K)*X(M,K) 20 CONTINUE CM(M)=CM(M)/WM 19 CONTINUE C C...CONVERT COORDINATES OF EACH ATOM TO CARTESIANS WITH ORIGIN AT C CENTRE OF MASS C DO 21 I=1,3 DO 21 J=1,NATPA X(I,J)=X(I,J)-CM(I) 21 CONTINUE C C...CONVERT TO PRINCIPAL AXES C 1. DIAGONALISE INERTIA TENSOR C (BORROWING V AND EIVR ARRAYS) C DO 80 I=1,3 DO 80 J=1,3 V(I,J)=0.0 EIVR(I,J)=0.0 80 CONTINUE DO 81 I=1,NAT V(1,1)=V(1,1)+WT(I)*(X(2,I)**2+X(3,I)**2) V(2,2)=V(2,2)+WT(I)*(X(1,I)**2+X(3,I)**2) V(3,3)=V(3,3)+WT(I)*(X(1,I)**2+X(2,I)**2) V(1,2)=V(1,2)-WT(I)*X(1,I)*X(2,I) V(1,3)=V(1,3)-WT(I)*X(1,I)*X(3,I) V(2,3)=V(2,3)-WT(I)*X(2,I)*X(3,I) 81 CONTINUE V(2,1)=V(1,2) V(3,1)=V(1,3) V(3,2)=V(2,3) C CALL DIAG(3,1) ----- DO 185 I=1,3 DO 185 J=1,3 ROTMAT(I,J)=EIVR(I,J) 185 CONTINUE C C...2.ROTATE CARTESIANS: X = ROTMAT (3x3) * X (3xNAT) C (note that the new X,Y,Z are inertial axes A,B,C in that order) C C...check sign of determinant to establish whether handedness of axis C system not changed, take remedial action if necessary rotdet= * rotmat(1,1)*(rotmat(2,2)*rotmat(3,3)-rotmat(3,2)*rotmat(2,3)) * -rotmat(1,2)*(rotmat(2,1)*rotmat(3,3)-rotmat(3,1)*rotmat(2,3)) * +rotmat(1,3)*(rotmat(2,1)*rotmat(3,2)-rotmat(3,1)*rotmat(2,2)) if(rotdet.lt.0.d0)then write(iw,'(1x/'' -> Sign of second row of rotation matrix'', * '' reversed to preserve handedness'')') rotmat(1,2)=-rotmat(1,2) rotmat(2,2)=-rotmat(2,2) rotmat(3,2)=-rotmat(3,2) endif C C...read external rotation matrix if specified - this is to be selected so C that handedness of coordinates is not changed, ie determinant is +1. C Note that for C3v molecules conventionally the C3 (z) axis is with the C threefold group at the negative end, the xz plane is a symmetry plane, C and the off-axis atom on that plane is for -ve x c IF(INPROT.EQ.1)THEN write(iw,'(1x/'' NOTE: External rotation matrix -->'')') DO 884 J=1,3 READ(IR,*)(ROTMAT(I,J),I=1,3) 884 CONTINUE ENDIF C C...rotate DO 82 I=1,NATPA XTEMP=0.0 YTEMP=0.0 ZTEMP=0.0 DO 83 J=1,3 XTEMP = XTEMP+ROTMAT(J,1)*X(J,I) YTEMP = YTEMP+ROTMAT(J,2)*X(J,I) ZTEMP = ZTEMP+ROTMAT(J,3)*X(J,I) 83 CONTINUE X(1,I)= XTEMP X(2,I)= YTEMP X(3,I)= ZTEMP 82 CONTINUE WRITE(IW,206) DO 86 J=1,3 WRITE(IW,207)LABEL(J),(ROTMAT(I,J),I=1,3) 86 CONTINUE rotdet= * rotmat(1,1)*(rotmat(2,2)*rotmat(3,3)-rotmat(3,2)*rotmat(2,3)) * -rotmat(1,2)*(rotmat(2,1)*rotmat(3,3)-rotmat(3,1)*rotmat(2,3)) * +rotmat(1,3)*(rotmat(2,1)*rotmat(3,2)-rotmat(3,1)*rotmat(2,2)) write(iw,227)rotdet 227 format(' determinant=',f9.6) 206 FORMAT(/' ROTATION MATRIX TO PRINCIPAL COORDINATES:'/ * ' -----------------------------------------' */16X,'XOLD',16X,'YOLD',16X,'ZOLD'/) 207 FORMAT(1X,A,3(F20.6)) 778 CONTINUE C C...3.RESTORE V AND EIVR C DO 84 I=1,3 DO 84 J=1,3 V(I,J)=0.0 EIVR(I,J)=0.0 84 CONTINUE WRITE(IW,205) DO 188 K=1,NATPA WRITE(IW,204)K,(X(KK,K),KK=1,3),WT(K) 188 CONTINUE 101 FORMAT(1X,A4,2I5,3F10.4,2I5,F10.4) 201 FORMAT(/' ATOMIC COORDINATE INPUT (Angstr):'/ * ' ---------------------------------'/ * 5X,'ATOM',2X,'ORIGIN',2X, *'POLAR',3X,'R',9X,'THETA',9X,'PHI',5X,'LL',4X,'NAX',4X,'ALFA'/ * 20X,'AXIS') 202 FORMAT(1X,I2,2X,A4,I5,4X,A1,F12.6,F12.6,F12.6,2I5,F12.6) 203 FORMAT(/' CARTESIAN COORDINATES (Angstr):'/ * ' -------------------------------' * /' ATOM NO.',9X,'X',13X,'Y',13X,'Z',9X,'MASS (u)'/) 205 FORMAT(/' PRINCIPAL COORDINATES (Angstr):'/ * ' ------------------------------' * /' ATOM NO.',9X,'a',13X,'b',13X,'c',9X,'MASS (u)'/) 204 FORMAT(1X,I2,7X,3(F12.6,2X),F12.6) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE VIBFR(NAT,NFK,IFVCT,NMOD,IPRFUL,fscale) C C VIBFR SETS UP AND SOLVES THE VIBRATIONAL SECULAR EQUATIONS C C NMOD - the number of normal modes (and non-redundant symmetry coordinates) C NINTC - the number of internal coordinates C NAT - the number of atoms C N - the number of Cartesian coordinates = 3*NAT C X - Cartesian coordinates (3 x NAT) C C V - (NxN) array in which the G matrix is set up in Cartesian C coordinates and then diagonalised. C After diagonalisation the eigenvalues are in V(i,i) in C increasing order of magnitude and i=1,6 are six rotations C and translations C EIVR - eigenvectors from the diagonalisation, inclusive of C rotations and translations C EU - unweighted eigenvectors (inclusive of tr+rot) C EM - mass weighted eigenvectors (inclusive of tr+rot) C EL - just the normal mode mass-weighted eigenvecotrs (NxNMOD) C C GM - used initially for the G matrix and then for the L matrix C BM - B matrix C G - vector of square roots of masses of the atoms (3*NAT) C FINTC - internal force constant matrix (NINTCxNINTC) for PED calc. C RLAMBD(i,1) - eigenvalues transferred from V(i,i) C RLAMBD(i,2) - vibrational frequencies (cm-1) C INTERN(i,j) - definitions of internal coordinates, j=1,4 atoms, C j=5 coordinate type, j=6 unique index C NATOMS(i) - the number of defining atoms for coordinate type i C COORD(i) - descriptive string for coordinate type i C C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/V(maxmod,maxmod),EIVR,WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) c REAL*8 EIVR(maxmod,maxmod),EL(maxmod,maxmod), * EM(maxat,3),EU(maxat,3) c REAL*4 BM(maxmod,maxmod),GM(maxmod,maxmod), * SCALL(maxmod),U(maxmod,maxmod), * SYML(maxmod,maxmod),FINTC(maxmod,maxmod) REAL*4 SIGMA(maxmod,maxmod),LPRIME(maxmod,maxmod), * delta(maxmod,maxmod),work(maxmod,maxmod),work1(maxmod,maxmod), * del298(maxmod,maxmod),SIG298(maxmod,maxmod),kappal c CHARACTER COORD(7)*4 c DIMENSION B(72),DC(12),C1(3),C2(3),C3(3),LIST(12,2),RR(3), * IBAS(6,maxat),MULT(maxat),G(maxmod), * NATOMS(6),INTERN(maxmod,6),listb(4,2),indexb(maxat) INTEGER*2 IEQUIV(maxmod) CHARACTER HEADER*15,TEXT(6) DATA TEXT/'i',',','i','i',',','j'/ DATA COORD/' ','STR.','BEND','OPLA','LIBE','TORS','tors'/ DATA ZERO/0.0D0/,ONE/1.0D0/,CONV/1302.7811D0/ DATA NATOMS/2,3,4,4,2,4/ NINTC=0 C C...G - THE ARRAY OF SQUARE ROOTS OF ATOMIC MASSES C DO 1 I=1,NAT DO 2 IB=1,6 IBAS(IB,I)=ZERO 2 CONTINUE MULT(I)=0 J=3*I-2 G(J)=DSQRT(WT(I)) G(J+1)=G(J) G(J+2)=G(J) 1 CONTINUE C C...V - THE POTENTIAL ENERGY MATRIX C N=3*NAT DO 3 I=1,N DO 3 J=1,N V(I,J)=ZERO BM(I,J)=ZERO fintc(i,j)=0.d0 3 CONTINUE C write(iw,2201)fscale 2201 format(1x//' NOTE: force constants are scaled by factor of', * f6.3,' relative to input.') WRITE(IW,201) C C...FK - THE ARRAY OF FORCE CONSTANT VALUES, READ THESE IN 10 PER LINE C in units of mdyn/A, mdyn/rad, mdyn A/rad^2 for stretch, stretch-bend, C bend resp. C NTIMES=NFK/10 IF(NFK-NTIMES*10.GT.0)NTIMES=NTIMES+1 L=1 DO 499 I=1,NTIMES J=L+9 IF(J.GT.NFK)J=NFK READ(IR,*) (FK(INDEX),INDEX=L,J) L=L+10 499 CONTINUE if(fscale.ne.1.0D0)then do 5501 i=1,nfk fk(i)=fscale*fk(i) 5501 continue endif C C C...READ A CARD SPECIFYING TWO INTERNAL COORDINATES CONTRIBUTING A C TERM TO THE POTENTIAL ENERGY: C 4 READ(IR,102)ICC1,(LIST(I,1),I=1,4),INDEX,ICC2,(LIST(I,2),I=1,4), * HEADER C C...CHECK FOR FINAL DATA CARD BEFORE CALCULATING FREQUENCIES C IF(ICC1.EQ.0) GO TO 18 C C...FISH OUT DEFINITIONS OF INTERNAL COORDINATES C (this is only required for PED calculation and printout of various C matrices) C NOTE: sixth column of INTERN contains INDXI which uniquely identifies C the internal mode, and can be used for assigning off-diagonal force C constants. C IF(ICC2.EQ.0)THEN NINTC=NINTC+1 NUM=NATOMS(ICC1) DO 31 J=1,4 INTERN(NINTC,J)=LIST(J,1) 31 CONTINUE INTERN(NINTC,5)=ICC1 INDXI=LIST(1,1) DO 32 J=2,NUM INDXI=INDXI+MAXAT**(J-1)*LIST(J,1) 32 CONTINUE INTERN(NINTC,6)=INDXI ELSE do 757 j=1,4 listb(j,1)=list(j,1) listb(j,2)=list(j,2) 757 continue ENDIF INTFC=0 C C...CHECK WHETHER CONTRIBUTION TO POTENTIAL ENERGY IS A SQUARE TERM C OR A CROSS TERM : C INTFC = 0 IF FK(INDEX) IS A DIAGONAL FORCE CONSTANT C INTFC = 1 IF FK(INDEX) IS AN OFF-DIAGONAL (INTERACTION) F.C. C IF(ICC2.NE.0) INTFC=1 KC1=ICC1+1 KC2=ICC2+1 KW=1 IF(INTFC.EQ.1) KW=4 WRITE(IW,202) TEXT(KW),TEXT(KW+1),TEXT(KW+2),INDEX,FK(INDEX), * COORD(KC1),(LIST(I,1),I=1,4),COORD(KC2),(LIST(I,2),I=1,4), * HEADER M=0 NORM=1 C C...CONVERT THE TWO INTERNAL COORDINATES INVOLVED IN THE POTENTIAL TERM TO C CARTESIAN DISPLACEMENT COORDINATES FOR EACH ATOM USING ANALYTICAL C EXPRESSIONS FOR S VECTORS. C (ELEMENTS OF ARRAY B ARE COMPONENTS OF S VECTORS) C DO 5 I=1,2 IT=I ITYP=ICC1 IF(I.EQ.2.AND.INTFC.EQ.0) GO TO 7 IF(I.EQ.2) ITYP=ICC2 GO TO (100,200,300,400,500,600), ITYP C C...STRETCHING COORDINATE FOR BOND P-Q SPECIFIED BY : C LIST(1,I) = P , LIST(2,I) = Q C 100 NUM=2 CALL DCS(IT,DC,RR,1,LIST,ITYP) ----- DO 110 L=1,3 M=M+1 B(M)=DC(L) B(M+3)=-DC(L) 110 CONTINUE C C...CONSTRUCT A TOPOLOGICAL MAP OF THE BONDED ATOMS C C IBAS is a (6xNATOMS) array in which information concerning connectivity C of the molecule is store ie. numbers of atoms to which a given C atom is bonded C MULT is a NATOMS length vector containing the multiplicity of each atom C ie. the number of atoms each atom is bonded to C IA=LIST(1,I) IB=LIST(2,I) IF(MULT(IA).EQ.0) GO TO 120 DO 130 LA=1,NLA IF(IBAS(LA,IA).EQ.IB) GO TO 6 130 CONTINUE 120 MULT(IA)=MULT(IA)+1 MULT(IB)=MULT(IB)+1 NLA=MULT(IA) NLB=MULT(IB) IBAS(NLA,IA)=IB IBAS(NLB,IB)=IA GO TO 6 C C...BENDING COORDINATE FOR ANGLE P-Q-R SPECIFIED BY : C LIST(1,I) = P , LIST(2,I) = Q , LIST(3,I) = R C 200 NUM=3 NT=LIST(3,I) LIST(3,I)=LIST(2,I) LIST(2,I)=NT CALL DCS(IT,DC,RR,2,LIST,ITYP) ----- CIJK=DOT(1,2,DC,DC) CALL CROSS(1,2,DC,DC,C1,SIJK) ----- IF(SIJK-0.0001D0) 700,700,210 210 DO 220 J=1,3 M=M+1 B(M)=(CIJK*DC(J)-DC(J+3))/(RR(1)*SIJK) B(M+3)=(CIJK*DC(J+3)-DC(J))/(RR(2)*SIJK) B(M+6)=-(B(M)+B(M+3)) 220 CONTINUE GO TO 6 C C...OUT-OF-PLANE BENDING COORDINATE FOR FOUR COPLANAR ATOMS, C APICAL ATOM Q BONDED TO END ATOM P AND ANCHOR ATOMS R1 AND R2, C SPECIFIED BY: LIST(1,I)=P ,LIST(2,I)=R1 ,LIST(3,I)=R2 ,LIST(4,I)=Q C ie. P is attached to Q and is bending out of the plane of R1-Q-R2 C NOTE: it appears that for two OPLA coordinates with central atoms C at opposite ends of a bond, R1 and R2 in both coordinates have to C be specified in the same sense - ie both pairs in either clockwise or C anticlockwise order. C 300 NUM=4 CALL DCS(IT,DC,RR,3,LIST,ITYP) ----- CALL CROSS(1,2,DC,DC,C3,R3) ----- CALL CROSS(3,1,DC,DC,C2,R2) ----- CALL CROSS(2,3,DC,DC,C1,R1) ----- DTKL=DOT(2,3,DC,DC) DO 310 J=1,3 DC(9+J)=C1(J) 310 CONTINUE DTI1=DOT(1,4,DC,DC) IF(R1-0.0001D0) 700,700,320 320 SINT=DTI1/R1 COST=DSQRT(ONE-SINT*SINT) IF(COST-0.0001D0) 700,700,330 330 TANT=SINT/COST CTR1=COST*R1 DO 340 J=1,3 M=M+1 B(M)=(C1(J)*R1/CTR1 - TANT*DC(J))/RR(1) B(M+3)=(C2(J)*R2/CTR1 1 - TANT*(DC(J+3)-DTKL*DC(J+6))/(R1*R1))/RR(2) B(M+6)=(C3(J)*R3/CTR1 1 - TANT*(DC(J+6)-DTKL*DC(J+3))/(R1*R1))/RR(3) B(M+9)=-(B(M)+B(M+3)+B(M+6)) 340 CONTINUE GO TO 6 C C...LINEAR BENDING COORDINATE FOR COLLINEAR ATOMS P-Q-R IN PLANE OF C DUMMY ATOM D SPECIFIED BY : C LIST(1,I) = P , LIST(2,I) = Q , LIST(3,I) = R , LIST(4,I) = D C 400 NUM=3 CALL DCS(IT,DC,RR,3,LIST,ITYP) ----- RIIJ=ONE/RR(1) RIJK=ONE/RR(2) DO 410 J=1,3 M=M+1 B(M)=-DC(J+6)*RIIJ B(M+6)=-DC(J+6)*RIJK B(M+3)=-(B(M)+B(M+6)) 410 CONTINUE GO TO 6 C C...GENERAL TORSIONAL INTERNAL COORDINATE ABOUT BOND J-K C 500 CALL TWIST(IT,M,LIST,B,NORM,NUM,IBAS,MULT) ----- GO TO 6 C C...SIMPLE TORSIONAL INTERNAL COORDINATE I-J-K-L C 600 NUM=4 CALL STWIST(IT,M,LIST,B) ----- C 6 M=36 IF(I.EQ.1) M1=NUM IF(I.EQ.2) M2=NUM GO TO 5 C C...B MATRIX ELEMENTS FOR SECOND COORDINATE OF A SQUARE TERM C (DIAGONAL FORCE CONST.) ARE EQUAL TO THOSE OF THE FIRST COORDINATE C 7 NUM3=NUM*3 DO 8 K=1,NUM3 B(36+K)=B(K) 8 CONTINUE M2=M1 DO 9 K=1,NUM LIST(K,2)=LIST(K,1) 9 CONTINUE IF(ITYP.EQ.5) NORM=NORM*NORM C 5 CONTINUE C C...Set up internal coordinate force constant matrix FINTC for PED C F=FK(INDEX)/NORM c if(icc2.eq.0)then fintc(NINTC,NINTC)=F else indx1=listb(1,1) indx2=listb(1,2) do 752 j=2,4 INDX1=INDX1+MAXAT**(J-1)*LISTB(J,1) INDX2=INDX2+MAXAT**(J-1)*LISTB(J,2) 752 continue do 753 j=1,nintc if(indx1.eq.intern(j,6))then intf1=j goto 754 endif 753 continue write(*,'(a)')' !!!!!! Internal coord indexing problem !!!!!!' 754 do 755 j=1,nintc if(indx2.eq.intern(j,6))then intf2=j goto 756 endif 755 continue write(*,'(a)')' !!!!!! Internal coord indexing problem !!!!!!' 756 fintc(intf1,intf2)=F fintc(intf2,intf1)=F endif C C...CONSTRUCT POTENTIAL ENERGY MATRIX IN CARTESIAN COORDINATES : EACH C CONTRIBUTION IS A PRODUCT OF A FORCE CONSTANT WITH S VECTOR C COMPONENTS FROM THE TWO INTERNAL COORDINATES C 10 IA=0 DO 11 IK=1,M1 DO 11 IX=1,3 IA=IA+1 IM=3*(LIST(IK,1)-1)+IX C C...the following line sets up the full B matrix (for printout and PED C calculation only) IF(ICC2.EQ.0)BM(NINTC,IM)=B(IA) C DR=B(IA)*F IB=36 DO 11 JK=1,M2 DO 11 JX=1,3 IB=IB+1 JM=3*(LIST(JK,2)-1)+JX DELV=DR*B(IB) V(IM,JM)=V(IM,JM)+DELV IF(INTFC.EQ.1) V(JM,IM)=V(JM,IM)+DELV 11 CONTINUE GO TO 4 C C...CONVERT TO MASS-WEIGHTED CARTESIAN COORDINATES C 18 DO 12 I=1,N DO 12 J=1,N V(I,J)=V(I,J)/(G(I)*G(J)) 12 CONTINUE C C...OBTAIN EIGENVALUES OF SYMMETRIC MATRIX V BY DIAGONALIZATION USING C A SIMILARITY TRANSFORMATION WITH ORTHOGONAL EIGENVECTOR MATRIX C (eigenvalues are returned on the diagonal of V, eigenvectors (matrix l) C are returned columnwise in EIVR, so that rows refer to internal C coordinate and columns to normal coordinate - first five or six are C translations and rotations) C CALL DIAG(N,1) ----- C C...CONVERT EIGENVALUES TO FREQUENCIES IN CM**-1 C 984 DO 13 I=1,N TEMP=V(I,I) FREQ=DSQRT(DABS(TEMP))*CONV V(I,I)=DSIGN(FREQ,TEMP) 13 CONTINUE C C C...PRINT OUT NORMAL COORDINATES IN CARTESIANS IF IFVCT=1, OR ONLY NORMAL C MODE FREQUENCIES IF IFVCT=0 C EM and EU are eigenvectors for a given mode (NATOMSx3), EM is the mass C weighted eigenvector, EU is the unweighed eigenvector. C 14 IF(IFVCT.EQ.1)THEN WRITE(IW,204) ELSE WRITE(IW,108) ENDIF WRITE(*,203) (V(I,I),I=1,N) C NMOD=0 DO 16 I=1,N UMAG=ZERO DO 17 J=1,NAT NU=3*(J-1) GJ=G(NU+1) DO 17 K=1,3 EMJK=EIVR(NU+K,I) EM(J,K)=EMJK EUJK=EMJK/GJ EU(J,K)=EUJK UMAG=UMAG+EUJK*EUJK 17 CONTINUE UMAG=DSQRT(UMAG) DO 20 J=1,NAT DO 20 K=1,3 EU(J,K)=EU(J,K)/UMAG 20 CONTINUE C C...EL contains only the normal mode mass-weighted eigenvectors C IF(V(I,I).GT.1.D0)THEN NMOD=NMOD+1 RLAMBD(NMOD,1)=(V(I,I)/CONV)**2 RLAMBD(NMOD,2)=V(I,I) II=0 DO 38 JJ=1,NAT DO 38 JJJ=1,3 II=II+1 EL(II,NMOD)=EM(JJ,JJJ) 38 CONTINUE ENDIF C IF(IFVCT.EQ.1)THEN WRITE(IW,205) NMOD,V(I,I),(M,(EM(M,L),L=1,3), * (EU(M,L),L=1,3),M=1,NAT) ELSE WRITE(IW,109)NMOD,V(I,I) ENDIF 16 CONTINUE C C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C...VARIOUS OPTIONAL MATRICES: C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...B matrix C IF(IPRFUL.GT.0)THEN iprcol=15 iprat=iprcol/3 c WRITE(IW,514) 514 FORMAT(//1x,78(1h-)/' B MATRIX (NINTx3N):'/1x,78(1h-)/) c iprs=1 iprf=N if(iprf.gt.iprcol)iprf=iprcol natst=1 natfi=NAT if(natfi.gt.iprat)natfi=iprat c 515 WRITE(IW,44)(' a.',I,' b.',I,' c.',I,I=natst,natfi) 44 FORMAT(19x,5(A,I2,A,I2,A,I2)) WRITE(IW,'(1X)') c DO 45 I=1,NINTC WRITE(IW,46)COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (BM(I,J),J=iprs,iprf) 45 CONTINUE 46 FORMAT(1X,A4,I3,3(',',I2),' ',15F7.4) c if(iprf.lt.N)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.N)iprf=N natst=natst+iprat natfi=natfi+iprat if(natfi.gt.NAT)natfi=NAT WRITE(IW,'(1X)') goto 515 endif c ENDIF C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C -1 C...G matrix = B M B' C DO 66 I=1,NINTC DO 66 J=1,NINTC UMAG=ZERO DO 51 JJ=1,N UMAG=UMAG+BM(I,JJ)*BM(J,JJ)/(G(JJ))**2 51 CONTINUE GM(I,J)=UMAG 66 CONTINUE C IF(IPRFUL.GT.0)THEN if(iprful.eq.1)iprcol=17 if(iprful.eq.2)iprcol=10 c WRITE(IW,64) 64 FORMAT(//1x,78(1h-)/' G MATRIX = B M^-1 B'' (NINTxNINT):' * /1x,78(1h-)) c iprs=1 iprf=nintc if(iprf.gt.iprcol)iprf=iprcol c 510 if(iprful.eq.1)then WRITE(IW,501)(I,I=iprs,iprf) else WRITE(IW,502)(I,I=iprs,iprf) endif 501 FORMAT(2X,'Internal coordinate ',17i6) 502 FORMAT(2X,'Internal coordinate ',10i10) WRITE(IW,'(1x)') c DO 67 I=1,NINTC if(iprful.eq.1)then WRITE(IW,503)I,COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (GM(I,J),J=iprs,iprf) else WRITE(IW,504)I,COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (GM(I,J),J=iprs,iprf) endif 67 CONTINUE 503 FORMAT(i3,' = ',A4,I3,3(',',I2),' ',17F6.3) 504 FORMAT(i3,' = ',A4,I3,3(',',I2),' ',10F10.6) c if(iprf.lt.nintc)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.nintc)iprf=nintc WRITE(IW,'(1X)') goto 510 endif ENDIF C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C -1/2 C...L matrix = B M l (l = Lscript = the conversion matrix from normal to C mass weighted coordinates, as in q = l Q. No C normalization problems as to magnitude occur for l C since it is unity normalized from consideration C of kinetic energy, and it comes out that way from C the diagonalization routine) C DO 70 I=1,NINTC DO 70 J=1,NMOD UMAG=0.D0 DO 71 JJ=1,N UMAG=UMAG+BM(I,JJ)*EL(JJ,J)/G(JJ) 71 CONTINUE GM(I,J)=UMAG 70 CONTINUE C IF(IPRFUL.GT.0)THEN if(iprful.eq.1)iprcol=17 if(iprful.eq.2)iprcol=10 c WRITE(IW,65) 65 FORMAT(//1x,78(1h-)/ * ' L MATRIX = B M^-1/2 l (NINTxNMODES, normalized', * ' to G = LL''):'/1x,78(1h-)) c iprs=1 iprf=nmod if(iprf.gt.iprcol)iprf=iprcol c 505 if(iprful.eq.1)then WRITE(IW,506)(I,I=iprs,iprf) WRITE(IW,507)(idnint(RLAMBD(I,2)),I=iprs,iprf) else WRITE(IW,508)(I,I=iprs,iprf) WRITE(IW,509)(RLAMBD(I,2),I=iprs,iprf) endif 506 FORMAT(2X,'Normal mode -> ',17i6) 507 FORMAT(18x, 17i6) 508 FORMAT(2X,'Normal mode -> ',10i10) 509 FORMAT(18x, 10F10.2) WRITE(IW,'(1x)') c DO 63 I=1,NINTC IF(IPRFUL.EQ.1)then WRITE(IW,511)COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (GM(I,J),J=iprs,iprf) ELSE WRITE(IW,512)COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (GM(I,J),J=iprs,iprf) endif 63 CONTINUE 511 FORMAT(A4,I3,3(',',I2),' ',17F6.3) 512 FORMAT(A4,I3,3(',',I2),' ',10F10.6) c if(iprf.lt.nmod)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.nmod)iprf=nmod WRITE(IW,'(1X)') goto 505 endif ENDIF c C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C c...Vibrational amplitudes: Cyvin eq.7.11,7.21,7.66 c C Sigma = L * delta * L' (nintc x nmod) (nmod x nmod) (nomd x nintc) c TWOKT=2.0*298.0*0.6950356 IF(IPRFUL.GT.0)THEN do 800 i=1,nmod do 800 j=1,nmod if(i.eq.j)then delta(i,j)=16.857631D0/rlambd(i,2) xx=rlambd(i,2)/TWOKT del298(i,j)=(16.857631D0/rlambd(i,2))*(1.0/tanh(xx)) else delta(i,j)=0.0 del298(i,j)=0.0 endif 800 continue c c delta*L' c do 801 i=1,nintc do 801 j=1,nmod lprime(i,j)=gm(j,i) 801 continue c call mmult(delta,Lprime,work,nmod,nmod,nintc) <----- call mmult(del298,Lprime,work1,nmod,nmod,nintc) <----- c c L * (delta * L') c call mmult(GM,work,sigma,nintc,nmod,nintc) <----- call mmult(GM,work1,sig298,nintc,nmod,nintc) <----- c c WRITE(IW,803) 803 FORMAT(//1x,78(1h-)/ * ' Mean square (Angstr^2) and mean amplitudes (Angstr) ', * 'for bonds:'/ * ' SIGMA = L delta L'' (NINTxNINT)'/55x, * 'Cyvin: eq.7.11,7.21,7.66' * /1x,78(1h-)// * ' T=0K T=298K'/ * ' u^2 sqrt(u^2) u^2 sqrt(u^2)'/) c DO 802 I=1,NINTC if(intern(i,5).eq.1) * WRITE(IW,512)COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * sigma(I,I),sqrt(abs(sigma(i,i))), * sig298(I,I),sqrt(abs(sig298(i,i))) 802 CONTINUE 818 FORMAT(A4,I3,3(',',I2),' ',2F10.6,3x,2F10.6) c ENDIF c C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C c...Vibrational amplitudes using R.Stolevik,H.M.Seip,S.J.Cyvin, c Chem.Phys.Lett. 15 (1972) 263-265 c IF(IPRFUL.GT.0)THEN WRITE(IW,809) 809 FORMAT(//1x,78(1h-)/ * ' Mean square (Angstr^2) and mean amplitudes (Angstr) for all ', * 'interatomic'/' separations:'/52x, * 'Stolevik: u^2 = eq.14,15,10'/64x,'K = eq.20,21,10' * /1x,78(1h-)//43x,'T=0K T=298K'/ * 24x,'R u^2 sqrt(u^2) K',12x, * 'u^2 sqrt(u^2) K'/) c c...go through the stretching coordinates c c dR d2R d2R c derivatives: -- = x/R --- = (R^2-x^2)/R^3 --- = - xy/R^3 c dx dx2 dxdy c c see Gwinn, eq.5 for first derivatives, eq.15 for second derivatives c nbonds=0 do 805 ii=1,nintc if(intern(II,5).ne.1)goto 805 na=intern(ii,1) nb=intern(ii,2) c nbonds=nbonds+1 if(na.lt.nb)then indexb(nbonds)=na*1000+nb else indexb(nbonds)=nb*1000+na endif c R=dsqrt( (x(1,na)-x(1,nb))**2 + (x(2,na)-x(2,nb))**2 + * (x(3,na)-x(3,nb))**2 ) c c...loop for u^2, Stolevik equations 14,15,10 c usq=0.d0 usq298=0.d0 do 806 l=1,nmod gammal=0.d0 do 808 i=1,nat if(i.ne.na.and.i.ne.nb)goto 808 ix=3*(i-1) do 810 ic=1,3 ix=ix+1 der=(x(ic,na)-x(ic,nb))/R if(i.eq.nb)der=-der ail=EL(ix,l)/G(ix) gammal=gammal+der*ail 810 continue 808 continue c usq=usq+gammal**2*delta(l,l) usq298=usq298+gammal**2*del298(l,l) 806 continue c c...loop for K, Stolevik equations 20,21,10 c rk=0.d0 rk298=0.d0 c do 822 l=1,nmod kappal=0.d0 c do 820 i=1,nat if(i.ne.na.and.i.ne.nb)goto 820 ix=3*(i-1) do 821 ic=1,3 ix=ix+1 do 824 j=1,nat if(j.ne.na.and.j.ne.nb)goto 824 jx=3*(j-1) do 825 jc=1,3 jx=jx+1 c if(ic.eq.jc)then x,x c c__x,x etc. der=(R**2-(x(ic,na)-x(ic,nb))**2)/R**3 x_na,x_na x_nb,y_nb if(ix.ne.jx)der=-der x_na,x_nb else c c__x,y etc. der=-(x(ic,na)-x(ic,nb))* * (x(jc,na)-x(jc,nb))/R**3 x_na,y_na x_nb,y_nb if(i.ne.j)der=-der x_na,y_nb x_nb,y_na endif c ail=EL(ix,l)/G(ix) ajl=EL(jx,l)/G(jx) kappal=kappal+0.5*der*ail*ajl 825 continue 824 continue 821 continue 820 continue c rk=rk+kappal*delta(l,l) rk298=rk298+kappal*del298(l,l) 822 continue c c WRITE(IW,819)COORD(INTERN(II,5)+1),(INTERN(II,J),J=1,4), * R,usq,sqrt(usq),RK, * usq298,sqrt(usq298),rk298 805 continue 819 FORMAT(A4,I3,3(',',I2),' ',F10.6,3x,3F10.6,3x,3F10.6) c c...go through nonbonded distances c write(iw,'(1x)') do 816 na=1,nat-1 do 811 nb=na+1,nat ii=na*1000+nb do 817 jj=1,nbonds if(ii.eq.indexb(jj))goto 811 817 continue c R=dsqrt( (x(1,na)-x(1,nb))**2 + (x(2,na)-x(2,nb))**2 + * (x(3,na)-x(3,nb))**2 ) c c...loop for u^2 c usq=0.d0 usq298=0.d0 do 812 l=1,nmod gammal=0.d0 do 813 i=1,nat if(i.ne.na.and.i.ne.nb)goto 813 ix=3*(i-1) do 814 ic=1,3 ix=ix+1 der=(x(ic,na)-x(ic,nb))/R if(i.eq.nb)der=-der ail=EL(ix,l)/G(ix) gammal=gammal+der*ail 814 continue 813 continue c usq=usq+gammal**2*delta(l,l) usq298=usq298+gammal**2*del298(l,l) 812 continue c c...loop for K c rk=0.d0 rk298=0.d0 c do 830 l=1,nmod kappal=0.d0 c do 829 i=1,nat if(i.ne.na.and.i.ne.nb)goto 829 ix=3*(i-1) do 828 ic=1,3 ix=ix+1 do 827 j=1,nat if(j.ne.na.and.j.ne.nb)goto 827 jx=3*(j-1) do 826 jc=1,3 jx=jx+1 c if(ic.eq.jc)then x,x c c__x,x etc. der=(R**2-(x(ic,na)-x(ic,nb))**2)/R**3 x_na,x_na x_nb,y_nb if(ix.ne.jx)der=-der x_na,x_nb else c c__x,y etc. der=-(x(ic,na)-x(ic,nb))* * (x(jc,na)-x(jc,nb))/R**3 x_na,y_na x_nb,y_nb if(i.ne.j)der=-der x_na,y_nb x_nb,y_na endif c ail=EL(ix,l)/G(ix) ajl=EL(jx,l)/G(jx) kappal=kappal+0.5*der*ail*ajl 826 continue 827 continue 828 continue 829 continue c rk=rk+kappal*delta(l,l) rk298=rk298+kappal*del298(l,l) 830 continue c c WRITE(IW,815)na,nb,R,usq,sqrt(usq),rk, * usq298,sqrt(usq298),RK298 811 continue 816 continue 815 format(4x,i3,' ...',i2,5x,F10.6,3x,3F10.6,3x,3F10.6) c ENDIF C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SYMMETRY COORDINATES C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c read(ir,'(a)')header if(header(1:8).ne.'U MATRIX')then backspace(ir) goto 957 endif DO 996 I=1,NMOD READ(IR,*)(U(I,J),J=1,NINTC) 996 CONTINUE C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c C... L.sym = U L C CALL MMULT(U,GM,SYML,NMOD,NINTC,NMOD) C c...normalization of eigenvectors so that elements of L.sym with largest C absolute value are made positive - ie for positive displacement of C dominant symmetry coordinate for each normal coordinate. C This is in line with the convention in Cyvin's book, p111, to ensure c reproducible relative phases of eigenvectors so that relative signs c of Coriolis constants have some meaning. The normalising coefficients c +-1 worked out for symmetric coordinate eigenvectors L.sym are then c applied to mass-weighted eigenvectors l. c C NOTE: 1/ note that all arrays except EIVR have rotation and translation C coordinates removed C 2/ in this version the number of symmetry coordinates is to be C the same as no of normal modes - ie redundant coordinates C not included C isumnm=0 isumsc=0 do 991 i=1,nmod vmax=0.0 do 992 ii=1,nmod if(abs(syml(ii,i)).gt.abs(vmax))THEN vmax=syml(ii,i) iequiv(i)=ii endif 992 continue isumnm=isumnm+i isumsc=isumsc+iequiv(i) scall(i)=1.0 if(vmax.lt.0.d0)then scall(i)=-1 do 993 ii=1,nmod syml(ii,i)=-syml(ii,i) 993 continue iii=n-nmod+i do 994 ii=1,n eivr(ii,iii)=-eivr(ii,iii) 994 continue endif 991 continue C C IF(IPRFUL.NE.0)THEN iprcol=17 C WRITE(IW,940) 940 FORMAT(//1x,78(1h-)/ * ' U matrix (NSYMxNINT, internals are ', * 'in the same order as in the B matrix):'/1x,78(1h-)/) c iprs=1 iprf=nintc if(iprf.gt.iprcol)iprf=iprcol c 517 write(IW,518)(I,I=iprs,iprf) 518 FORMAT(' internal->',17i6) write(IW,519) 519 FORMAT(' _symmetry'/'|') c DO 941 I=1,NMOD WRITE(IW,942)I,(U(I,J),J=iprs,iprf) 941 CONTINUE 942 FORMAT(I5,7x,17F6.3) c if(iprf.lt.nintc)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.nintc)iprf=nintc WRITE(IW,'(1X)') goto 517 endif C C- - - - C WRITE(IW,965) 965 FORMAT(//1x,78(1h-)/ * ' L.sym = U L.int (NSYMxNMOD, phases adjusted to make', * ' major elements +ve):'/1x,78(1h-)/) c iprs=1 iprf=nmod if(iprf.gt.iprcol)iprf=iprcol c 520 write(IW,518)(I,I=iprs,iprf) write(IW,519) c DO 963 I=1,NMOD WRITE(IW,942)I,(SYML(I,J),J=iprs,iprf) 963 CONTINUE write(iw,944)(int(scall(j)),j=iprs,iprf) 944 format(' phasing:',17i6) C if(iprf.lt.nmod)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.nmod)iprf=nmod WRITE(IW,'(1X)') goto 520 endif C C- - - - C write(iw,971) 971 format(//1x,78(1h-)/ * ' Correspondence between normal modes and symmetry', * ' coordinates:'/1x,78(1h-)/) c iprcol=30 iprs=1 iprf=nmod if(iprf.gt.iprcol)iprf=iprcol c 522 write(iw,969)' Normal mode: ',(i,i=iprs,iprf) write(iw,969)'Symmetry coord.: ',(iequiv(i),i=iprs,iprf) 969 format(1x,a,30i3) c if(iprf.lt.nmod)then iprs=iprs+iprcol iprf=iprf+iprcol if(iprf.gt.nmod)iprf=nmod WRITE(IW,'(1X)') goto 522 endif c if(isumnm.ne.isumsc)write(iw,970) 970 format(1x/' **** WARNING: a symmetry coordinate assigned ', * 'more than once'/) ENDIF C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...POTENTIAL ENERGY DISTRIBUTION: C C PED = Lambda^-1 J F_int [Califano], eq.8.6.14 C C J is matrix of derivatives of Lambda wrt F with elements: C C J_i_kl = L_ki L_li [Califano], eq.4.7.11, factor of 2 for C k.ne.l is not relevant here C 957 WRITE(IW,41) 41 FORMAT(//1x,78(1h-)/ * ' PED MATRIX = LAMBDA^-1 J F_int (NINTxNMODES):'// * ' (J is the Jacobian of LAMBDA wrt F_int)', * ' Califano: 8.6.14,4.7.11'/1x,78(1h-)) c DO 39 I=1,NMOD DO 40 J=1,NINTC V(I,J)=0.d0 DO 750 jj=1,nintc V(I,J)=V(I,J)+100.*FINTC(JJ,J)*GM(J,I)*GM(JJ,I)/RLAMBD(I,1) 750 CONTINUE 40 CONTINUE 39 CONTINUE c iprs=1 iprf=nmod if(iprf.gt.27)iprf=27 959 WRITE(IW,958)(I,I=iprs,iprf) 958 FORMAT(19X,50I4) WRITE(IW,141)(RLAMBD(I,2),I=iprs,iprf,2) WRITE(IW,142)(RLAMBD(I,2),I=iprs+1,iprf,2) 141 FORMAT(2X,'Normal mode -> ',20F8.1) 142 FORMAT(21X,20F8.1) WRITE(IW,'(1X)') c DO 43 I=1,NINTC WRITE(IW,42)COORD(INTERN(I,5)+1),(INTERN(I,J),J=1,4), * (IDNINT(V(J,I)),J=iprs,iprf) 43 CONTINUE if(iprf.lt.nmod)then iprs=iprs+27 iprf=iprf+27 if(iprf.gt.nmod)iprf=nmod WRITE(IW,'(1X)') goto 959 endif c 42 FORMAT(1X,A4,I3,3(',',I2),' ',50I4) C 101 FORMAT(7F10.6) 102 FORMAT(11I5,A15) 108 FORMAT(//1x,78(1h-)/' FREQUENCIES OF NORMAL MODES (cm-1):'/ * 1x,78(1h-)) 109 FORMAT(' Mode ',I2,':'F15.4) 201 FORMAT(//1x,78(1h-)/' SPECIFICATION OF MOLECULAR FORCE FIELD ', * '(in mdyn/A, mdyn/rad, mdyn A/rad^2):'/1x,78(1h-) * /4X,'FORCE CONSTANTS',24X,'COORDINATES' * /' TYPE INDEX VALUE',16X,'ATOMS',14X,'ATOMS'/) 202 FORMAT(1X,3A1,I5,3X,F10.6,6X,A4,4I3,3X,A4,4I3,A15) 203 FORMAT(/' FREQUENCIES OF NORMAL MODES (cm-1):'/ * ' -----------------------------------'/(6F12.4)) 204 FORMAT(//1x,78(1h-)/ *' FREQUENCIES (cm-1) AND CARTESIAN COORDINATES OF ', * 'NORMAL MODES'/1x,78(1h-)/ * /' FREQUENCY',10X,'MASS-WEIGHTED CARTESIANS',10X, * 'UNWEIGHTED CARTESIANS'/11X,'ATOM',6X,'a', 9X,'b', 9X,'c', * 12X,'a', 9X,'b',9X,'c'/) 205 FORMAT(1X,'Mode ',I2,':'/1X,F9.4,I4,2X,3F10.6,3X,3F10.6/ * (12X,I2,2X,3F10.6, 3X,3F10.6)) C GOTO 15 700 WRITE(IW,206) ITYP 206 FORMAT(' **** PROGRAM FAILS FOR COORDINATE TYPE',I3) C 15 RETURN END C C_____________________________________________________________________________ C SUBROUTINE TWIST(I,M,LIST,B,NORM,NUM,IBAS,MULT) C C TWIST CONVERTS A GENERAL TORSIONAL INTERNAL COORDINATE TO C CARTESIAN DISPLACEMENT COORDINATES FOR EACH ATOM C C I = whether coordinate is i or j (1 or 2) C LIST = atoms defining the coordinate C B = B matrix C MULT = no of atoms a given atom is bonded to C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) DIMENSION B(72),DC(12),C1(3),C2(3),RR(6), * LIST(12,2),LINK(12,2),MULT(maxat),IBAS(6,maxat) C DO 500 L=5,12 LIST(L,I)=0 500 CONTINUE DO 570 L=1,36 B(M+L)=0.0D0 570 CONTINUE C C...IA,IB ARE THE 'CENTRAL' ATOMS DEFINING THE BOND ABOUT WHICH C TORSION OCCURS C IA=LIST(1,I) IB=LIST(2,I) ICH=0 IN=3 MULTA=MULT(IA) MULTB=MULT(IB) C C...SELECT TWO 'TERMINAL' ATOMS C DO 510 LA=1,MULTA DO 510 LB=1,MULTB LINK(2,I)=IA LINK(3,I)=IB LINK(1,I)=IBAS(LA,IA) LINK(4,I)=IBAS(LB,IB) C C...CHECK THAT TERMINAL ATOMS ARE NOT CENTRAL ATOMS C IF(LINK(1,I).EQ.IB.OR.LINK(4,I).EQ.IA) GO TO 510 C C...CALC. DIRECTION COSINES,VECTOR AND SCALAR PRODUCTS OF BOND VECTORS C CALL DCS(I,DC,RR,3,LINK,5) ----- CALL CROSS(1,2,DC,DC,C1,SIJK) ----- CALL CROSS(3,2,DC,DC,C2,SJKL) ----- CIJK= -DOT(1,2,DC,DC) CJKL= -DOT(2,3,DC,DC) C C...CHECK FOR THE PRESENCE OF A COLLINEAR SUBSECTION : C IF SO THEN REJECT THIS CHAIN C IF(SIJK*SIJK-1.0D-04) 510,510,530 530 IF(SJKL*SJKL-1.0D-04) 510,510,520 520 ICH=ICH+1 C C...ICH - THE NUMBER OF FOUR-ATOM CHAINS NOT CONTAINING C COLLINEAR SUBSECTIONS C LINK(2,I)=LINK(4,I) C C...TEST WHETHER EACH TERMINAL ATOM HAS OCCURRED IN A PREVIOUS CHAIN : C IF NOT THEN ADD THE ATOM NUMBER TO LIST C DO 540 K=1,2 IZ=LINK(K,I) DO 550 L=3,IN LINK(K+2,I)=L IF(IZ.EQ.LIST(L,I)) GO TO 540 550 CONTINUE LIST(IN,I)=IZ IN=IN+1 540 CONTINUE C C...ITA,ITB - SPECIFY THE POSITION IN LIST OF THE TWO TERMINAL C ATOMS FOR THIS CHAIN C ITA=LINK(3,I) ITB=LINK(4,I) C C...CALCULATE COMPONENTS OF S VECTOR FOR EACH ATOM OF CHAIN C ADD CONTRIBUTIONS FROM THIS CHAIN TO B MATRIX ELEMENTS OBTAINED C FROM PREVIOUS CHAINS C DO 560 J=1,3 MJ=M+J DELS1= -C1(J)/(RR(1)*SIJK) DELS2= (C1(J)*(RR(2)-RR(1)*CIJK))/(RR(2)*RR(1)*SIJK) + * (CJKL*C2(J))/(RR(2)*SJKL) DELS4= -C2(J)/(RR(3)*SJKL) DELS3= -(DELS1+DELS2+DELS4) B(MJ)=B(MJ)+DELS2 B(MJ+3)=B(MJ+3)+DELS3 B(3*(ITA-1)+MJ)=B(3*(ITA-1)+MJ)+DELS1 560 B(3*(ITB-1)+MJ)=B(3*(ITB-1)+MJ)+DELS4 510 CONTINUE C C...NORM - NORMALISING COEFFICIENT C NUM - THE TOTAL NUMBER OF ATOMS IN LIST FOR THIS COORDINATE C NORM=NORM*ICH NUM=IN-1 C RETURN END C C_____________________________________________________________________________ C SUBROUTINE STWIST(I,M,LIST,B) C C STWIST CONVERTS THE SIMPLE SINGLE CHAIN TORSIONAL INTERNAL COORDINATE TO C CARTESIAN DISPLACEMENT COORDINATES FOR CONSTITUENT ATOMS C (adapted from TWIST) C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) DIMENSION B(72),DC(12),C1(3),C2(3),RR(6), * LIST(12,2),LINK(12,2) C c...conversion from a-b-c-d to b-c a d declaration consistent with logic of c TWIST nterm1=list(1,i) ncen1=list(2,i) ncen2=list(3,i) list(1,i)=ncen1 list(2,i)=ncen2 list(3,i)=nterm1 c DO 570 L=1,36 B(M+L)=0.0D0 570 CONTINUE C C...IA,IB ARE THE 'CENTRAL' ATOMS DEFINING THE BOND ABOUT WHICH C TORSION OCCURS C IA=LIST(1,I) IB=LIST(2,I) LINK(2,I)=IA LINK(3,I)=IB C C...THE TWO 'TERMINAL' ATOMS C LINK(1,I)=LIST(3,I) LINK(4,I)=LIST(4,I) C C...CALC. DIRECTION COSINES,VECTOR AND SCALAR PRODUCTS OF BOND VECTORS C CALL DCS(I,DC,RR,3,LINK,5) ----- CALL CROSS(1,2,DC,DC,C1,SIJK) ----- CALL CROSS(3,2,DC,DC,C2,SJKL) ----- CIJK= -DOT(1,2,DC,DC) ----- CJKL= -DOT(2,3,DC,DC) ----- C LINK(2,I)=LINK(4,I) C C...ITA,ITB - SPECIFY THE POSITION IN 'LIST' OF THE TWO TERMINAL C ATOMS FOR THIS CHAIN C ITA=3 ITB=4 C C...CALCULATE COMPONENTS OF S VECTOR FOR EACH ATOM OF CHAIN C ADD CONTRIBUTIONS FROM THIS CHAIN TO B MATRIX ELEMENTS OBTAINED C FROM PREVIOUS CHAINS C DO 560 J=1,3 MJ=M+J DELS1= -C1(J)/(RR(1)*SIJK) DELS2= (C1(J)*(RR(2)-RR(1)*CIJK))/(RR(2)*RR(1)*SIJK) + * (CJKL*C2(J))/(RR(2)*SJKL) DELS4= -C2(J)/(RR(3)*SJKL) DELS3= -(DELS1+DELS2+DELS4) B(MJ)=B(MJ)+DELS2 B(MJ+3)=B(MJ+3)+DELS3 B(3*(ITA-1)+MJ)=B(3*(ITA-1)+MJ)+DELS1 B(3*(ITB-1)+MJ)=B(3*(ITB-1)+MJ)+DELS4 560 CONTINUE C C RETURN END C C_____________________________________________________________________________ C SUBROUTINE DIAG(N,IEGEN) C C DIAG PERFORMS THE DIAGONALIZATION OF MATRIX A USING THE C VARIABLE THRESHOLD JACOBI METHOD C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/A(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) DATA ZERO/0.0D0/,ONE/1.0D0/ C 1 IF(IEGEN) 99,102,99 99 DO 101 J=1,N DO 100 I=1,N 100 EIVR(I,J)=ZERO 101 EIVR(J,J)=ONE C C...FIND ELEMENT OF A WITH LARGEST ABSOLUTE VALUE C 102 ATOP = ZERO DO 111 I=1,N DO 111 J=I,N IF(ATOP-DABS(A(I,J))) 104,111,111 104 ATOP=DABS(A(I,J)) 111 CONTINUE IF(ATOP)109,109,113 109 RETURN C C...CALCULATE THE STOPPING CRITERION -- DSTOP C 113 AVGF = (N*(N-1))*0.55D0 D = ZERO DO 114 JJ=2,N DO 114 II=2,JJ S=A(II-1,JJ)/ATOP 114 D=S*S+D DSTOP=(1.D-06)*D C C...CALCULATE THE THRESHOLD, THRSH C THRSH =DSQRT(D/AVGF)*ATOP C C...START A SWEEP C 115 IFLAG=0 DO 130 JCOL=2,N JCOL1=JCOL-1 DO 130 IROW=1,JCOL1 AIJ=A(IROW,JCOL) C C...COMPARE THE OFF-DIAGONAL ELEMENT WITH THRSH C IF(DABS(AIJ)-THRSH) 130,130,117 117 AII=A(IROW,IROW) AJJ=A(JCOL,JCOL) S=AJJ-AII C C...CHECK WHETHER CHOSEN ROTATION LESS THAN THE ROUNDING ERROR : C IF SO THEN DO NOT ROTATE C IF(DABS(AIJ)-1.D-09*DABS(S)) 130,130,118 118 IFLAG=1 C C...IF ROTATION VERY CLOSE TO 45 DEGREES, SET SIN AND COS = SQRT(O.5) C IF(1.D-10*DABS(AIJ)-DABS(S)) 116,119,119 119 S = 0.707106781D0 C=S GO TO 120 C C...CALCULATION OF SIN AND COS FOR ROTATION NOT VERY CLOSE TO 45 DEG. C COS = C , SIN = S C 116 T=AIJ/S S = 0.25D0/DSQRT(0.25D0 + T*T) C = DSQRT(0.5D0 + S) S = 2.0D0*T*S/C C C...CALCULATION OF THE NEW ELEMENTS OF MATRIX A C 120 DO 121 I=1,IROW T=A(I,IROW) U=A(I,JCOL) A(I,IROW)=C*T-S*U 121 A(I,JCOL)=S*T+C*U I2=IROW+2 IF(I2-JCOL)127,127,123 127 CONTINUE DO 122 I=I2,JCOL T=A(I-1,JCOL) U=A(IROW,I-1) A(I-1,JCOL)=S*U+C*T 122 A(IROW,I-1)=C*U-S*T 123 A(JCOL,JCOL)=S*AIJ+C*AJJ A(IROW,IROW)=C*A(IROW,IROW)-S*(C*AIJ-S*AJJ) DO 124 J=JCOL,N T=A(IROW,J) U=A(JCOL,J) A(IROW,J)=C*T-S*U 124 A(JCOL,J)=S*T+C*U C C...ROTATION COMPLETED. CHECK WHETHER EIGENVECTORS REQUIRED C IF(IEGEN) 131,126,131 131 DO 125 I=1,N T=EIVR(I,IROW) EIVR(I,IROW)=C*T-EIVR(I,JCOL)*S 125 EIVR(I,JCOL)=S*T+EIVR(I,JCOL)*C C C...CALCULATE THE NEW NORM D AND COMPARE WITH DSTOP C 126 CONTINUE S=AIJ/ATOP D=D-S*S IF(D-DSTOP)1260,129,129 C C...RECALCULATE DSTOP AND THRSH TO DISCARD ROUNDING ERRORS C 1260 D=ZERO DO 128 JJ=2,N DO 128 II=2,JJ S=A(II-1,JJ)/ATOP 128 D=S*S+D DSTOP=(1.D-06)*D 129 THRSH=DSQRT(D/AVGF)*ATOP 130 CONTINUE IF(IFLAG)115,134,115 C C...ARRANGE EIGENVALUES IN ORDER OF INCREASING MAGNITUDE C AND ASSOCIATED EIGENVECTORS IN THE SAME ORDER C 134 NU=N AMIN = ZERO DO 11 I=1,N IF(I.GE.NU) RETURN AMIN=A(I,I) J=1 31 IF (J .LT. I) GO TO 10 IF(A(J,J).GE.AMIN) GO TO 10 C C...IF IEGEN = -1, EXCLUDE UNCONVERGED EIGENVALUES FROM THIS ORDERING C TE=DABS(EIVR(N,J))+DABS(EIVR(N-1,J)) IF ((TE .GT. 0.05D0) .AND. (IEGEN .EQ. -1)) GO TO 15 II=I AMIN =A(J,J) A(J,J)=A(I,I) A(I,I)=AMIN 16 DO 12 K=1,N TEMP=EIVR(K,II) EIVR(K,II)=EIVR(K,J) 12 EIVR(K,J)=TEMP GO TO 10 15 AM=A(J,J) A(J,J)=A(NU,NU) A(NU,NU)=AM II=NU NU=NU-1 GO TO 16 10 J = J + 1 IF (J .LE. NU) GO TO 31 11 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ROTCS(LL,NAX,ALFA,S) C C ROTCS SETS UP AN ORTHOGONAL TRANSFORMATION MATRIX CORRESPONDING C TO A ROTATION BY AN ANGLE ALFA DEGREES ABOUT AXIS NAX C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S(3,3),T(3,3),ST(3,3) DATA ZERO/0.0D0/,ONE/1.0D0/,RAD/57.295779513082326D0/ ALFA=ALFA/RAD M=NAX K=MOD(NAX,3)+1 L=MOD(NAX+1,3)+1 C C...CHECK WHETHER THIS IS A ROTATION FROM THE ORIGINAL AXES (LL=1) C OR FROM AXES GENERATED BY AN IMMEDIATELY PREVIOUS ROTATION (LL=2) C IF(LL-1)1,2,3 1 K=1 L=2 M=3 ALFA=ZERO 2 S(K,K)=DCOS(ALFA) S(L,L)=S(K,K) S(M,M)=ONE S(L,K)=DSIN(ALFA) S(K,L)=-S(L,K) S(L,M)=ZERO S(K,M)=ZERO S(M,L)=ZERO S(M,K)=ZERO RETURN 3 T(K,K)=DCOS(ALFA) T(L,L)=T(K,K) T(M,M)=ONE T(L,K)=DSIN(ALFA) T(K,L)=-T(L,K) T(L,M)=ZERO T(K,M)=ZERO T(M,L)=ZERO T(M,K)=ZERO DO 4 J=1,3 DO 4 K=1,3 ST(J,K)=ZERO DO 4 L=1,3 4 ST(J,K)=ST(J,K)+S(J,L)*T(L,K) DO 5 J=1,3 DO 5 K=1,3 5 S(J,K)=ST(J,K) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE DCS(I,DC,RR,NBO,LINK,ITYP) C C DCS CALCULATES COMPONENTS OF UNIT VECTORS EII ALONG EACH BOND C DEFINING A SPECIFIED INTERNAL COORDINATE C C I - defines whether the first or second internal coordinate in a C potential term is to be processed (=1 or =2) C DC - vector for output of components of e.ii's, in order 1x,1y,1z,2x,2y... C RR - vector returning bond lengths of up to three bonds involved in the C internal coordinate C NBO - to be set on input to the number of bonds involved in the i.c. C LINK - array containing the definition of the internal coordinate C ITYP - the type of i.c. (str=1, bend=2, etc.) C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) DIMENSION LINK(12,2),DC(12),C(3),RR(3) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) C DO 1 IBO=1,NBO IM=LINK(IBO,I) IN=LINK(IBO+1,I) IF(ITYP.EQ.2) IN=LINK(3,I) IF(ITYP.EQ.3) IN=LINK(4,I) IF(ITYP.EQ.4.AND.IBO.EQ.3) IM=LINK(2,I) C C...CALCULATE COMPONENTS OF BOND VECTOR C DO 2 JJ=1,3 C(JJ)=X(JJ,IM)-X(JJ,IN) 2 CONTINUE C C...CALCULATE BOND LENGTH C R=DSQRT(C(1)*C(1)+C(2)*C(2)+C(3)*C(3)) RR(IBO)=R C C...CALCULATE DIRECTION COSINES C DO 1 LL=1,3 IK=3*(IBO-1)+LL DC(IK)=C(LL)/R 1 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE CROSS(II,JJ,VA,VB,VC,R) C C CROSS CALCULATES THE VECTOR PRODUCT OF TWO UNIT VECTORS EII X EJJ C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION VA(12),VB(12),VC(3) C I=3*II-2 J=3*JJ-2 VC(1)=VA(I+1)*VB(J+2)-VA(I+2)*VB(J+1) VC(2)=VA(I+2)*VB(J)-VA(I)*VB(J+2) VC(3)=VA(I)*VB(J+1)-VA(I+1)*VB(J) R=DSQRT(VC(1)*VC(1)+VC(2)*VC(2)+VC(3)*VC(3)) C C...R - THE SINE OF THE ANGLE BETWEEN THE TWO VECTORS C IF(R.LE.0.01D0) RETURN VC(1)=VC(1)/R VC(2)=VC(2)/R VC(3)=VC(3)/R C C...VC - UNIT VECTOR PERPENDICULAR TO PLANE OF EII AND EJJ C RETURN END C C_____________________________________________________________________________ C FUNCTION DOT(II,JJ,DC,DR) C C DOT IS THE SCALAR PRODUCT OF TWO UNIT VECTORS EII AND EJJ C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DC(12),DR(12) C I=3*II-2 J=3*JJ-2 DOT=DC(I)*DR(J)+DC(I+1)*DR(J+1)+DC(I+2)*DR(J+2) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE MMULT(A,B,C,N,M,L) C C C MULT WILL MULTIPLY TWO GENERAL, DOUBLE DIMENSIONED, DOUBLE PRECISION C ARRAYS AA AND BB, TO GIVE THE RESULT CC C C CC ( N TIMES L ) = AA ( N TIMES M ) * BB ( M TIMES L ) C C - A IS A TWO DIMENSIONAL ARRAY (NMAX x NMAX) CONTAINING C AA IN THE TOP LEFT HAND CORNER C C - B IS A TWO DIMENSIONAL ARRAY (NMAX x NMAX) CONTAINING C BB IN THE TOP LEFT HAND CORNER C C - C IS A TWO DIMENSIONAL ARRAY (NMAX x NMAX) CONTAINING C THE PRODUCT OF AA AND BB ON OUTPUT C C parameter (maxat=100,nmax=3*maxat) REAL*4 A(NMAX,NMAX),B(NMAX,NMAX),C(NMAX,NMAX),PROD C DO 1 I=1,N DO 1 J=1,L PROD=0.0D0 DO 2 K=1,M PROD=PROD+A(I,K)*B(K,J) 2 CONTINUE C(I,J)=PROD 1 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE VRCON(NAT,NMOD,IPRFUL) C C CALCULATION OF VARIOUS VIBRATION-ROTATION CONSTANTS C C NAT - the number of atoms C NMOD - the number of normal modes C IMPLICIT REAL*8 (A-H,O-Z) parameter (maxat=100,maxmod=3*maxat,maxfk=10000) COMMON /JOINT/V(maxmod,maxmod),EIVR(maxmod,maxmod),WT(maxat), * X(3,maxat),FK(maxfk),IR,IW,RLAMBD(maxmod,2),ROTMAT(3,3) COMMON /ZET/ZETA,CAAAK REAL*8 ZETA(maxmod,maxmod,3),CAAAK(maxmod,3) DIMENSION AI(3),RC(3),A(6,maxmod),SWT(maxat),DEPSIL(maxmod,3), * ALPHA(maxmod,3) DIMENSION TEMP(6),TAU(maxmod,21),DC(maxmod,5) DIMENSION CDOBL(5),CDPRO(5),CDEXP(5),ERRCD(5) CHARACTER SGN,HEADER*10 c DATA CONV/1302.7811D0/,CONV1/0.17222125D0/ DATA CONV2/16.857631D0/,CONV3/505379.01D0/, * CONV4/0.0225784/,RK/16.863D0/ DATA TAAAA,TBBBB,TCCCC,TAABB,TAACC,TBBCC,TABAB,TACAC,TBCBC, * TBBBC,TABBB,TAAAC,TAAAB,TACCC,TBCCC/15*0.D0/ C C...AI(3) - MOMENTS OF INERTIA (amu Angstr**2) C RC(3) - ROTATIONAL CONSTANTS (cm-1) C AA,BB,CC - ROTATIONAL CONSTANTS (MHz) C RLAMBD(i,1) - EIGENVALUES (i is a mode counter - rot+transl subtracted) C units: mdyne/(amu Angstr) C RLAMBD(i,2) - FREQUENCIES (cm-1) C A(N,MODE) - DERIVATIVES dI/dQ (amu^1/2 Angstr), N=1-6: AA,BB,CC,AB,AC,BC C C CONVERSION FACTORS: C C Q(RED) = 0.172219.(OMEGA**1/2).Q C C RC(X) /cm-1 = 16.858039/I(X) (I / amu Angstr+2) C /MHz = 505379.01/I(X) C C Tau /cm-1 = 0.0225784 * Tau / amu-2 Angstr-5 mdyn-1 C C Omega /cm-1 = 1302.781 ( Lambda )**(1/2) for Lambda evaluated from C GF with masses in u, and force constants C in mdyn/A, mdyn /rad, mdyn A/rad^2 ie. units C of Lambda: mdyn/(amu A) C CONV4=CONV4*29979.2458D0 ZERO=0.0D0 NMAX=maxmod N=3*NAT MDFRST=N-NMOD+1 C C...Moments of inertia are by default calculated from cartesian coordinates C for the atoms, but these can optionally be overriden by externally C supplied values C DO 1 I=1,3 AI(I)=0.0 1 CONTINUE DO 2 I=1,NAT SWT(I)=DSQRT(WT(I)) AI(1)=AI(1) + WT(I)*(X(2,I)**2+X(3,I)**2) AI(2)=AI(2) + WT(I)*(X(1,I)**2+X(3,I)**2) AI(3)=AI(3) + WT(I)*(X(1,I)**2+X(2,I)**2) 2 CONTINUE C C C...INERTIAL DERIVATIVES C DO 3 I=1,6 DO 3 J=1,NMAX A(I,J)=ZERO 3 CONTINUE II=0 DO 4 I=MDFRST,N II=II+1 DO 5 J=1,NAT DO 99 K=1,3 NJ=3*(J-1)+K TEMP(K)= X(K,J)*EIVR(NJ,I) 99 CONTINUE K=3*(J-1) TEMP(4)=X(1,J)*EIVR(K+2,I)+X(2,J)*EIVR(K+1,I) TEMP(5)=X(1,J)*EIVR(K+3,I)+X(3,J)*EIVR(K+1,I) TEMP(6)=X(2,J)*EIVR(K+3,I)+X(3,J)*EIVR(K+2,I) A(1,II)=A(1,II)+2.0*SWT(J)*(TEMP(2)+TEMP(3)) A(2,II)=A(2,II)+2.0*SWT(J)*(TEMP(1)+TEMP(3)) A(3,II)=A(3,II)+2.0*SWT(J)*(TEMP(1)+TEMP(2)) A(4,II)=A(4,II)-SWT(J)*TEMP(4) A(5,II)=A(5,II)-SWT(J)*TEMP(5) A(6,II)=A(6,II)-SWT(J)*TEMP(6) 5 CONTINUE 4 CONTINUE C C C...OUTPUT C pa=-0.5d0*(Ai(1)-Ai(2)-Ai(3)) pb=-0.5d0*(Ai(2)-Ai(1)-Ai(3)) pc=-0.5d0*(Ai(3)-Ai(1)-Ai(2)) WRITE(IW,206) WRITE(IW,207)'Internal:',AI(1),AI(2),AI(3) WRITE(IW,207)' ',pa,pb,pc dv=-2.d0*pc WRITE(IW,205)'Ic-Ia-Ib=',dv 206 FORMAT(//1x,78(1h-)/ * ' MOMENTS OF INERTIA I.a, I.b, I.c (u Angstr^2):'/ * ' PLANAR MOMENTS P.a, P.b, P.c (u Angstr^2):'/ * 1x,78(1h-)) 207 FORMAT(1X,A9,5X,3F20.8) 205 FORMAT(1x/46x,A9,F20.8) C read(ir,'(a)')header if(header.ne.'MOMENTS OF')THEN BACKSPACE(ir) goto 101 endif READ(IR,*,ERR=101,END=101)RIA,RIB,RIC AI(1)=RIA AI(2)=RIB AI(3)=RIC write(IW,'(1x)') WRITE(IW,207)'External:',AI(1),AI(2),AI(3) pa=-0.5d0*(Ai(1)-Ai(2)-Ai(3)) pb=-0.5d0*(Ai(2)-Ai(1)-Ai(3)) pc=-0.5d0*(Ai(3)-Ai(1)-Ai(2)) WRITE(IW,207)' ',pa,pb,pc dv=-2.d0*pc WRITE(IW,205)'Ic-Ia-Ib=',dv C 101 IF(AI(1).EQ.0.0D0)AI(1)=conv3*1.D-10 if(AI(2).eq.AI(3))AI(3)=1.0000001d0*AI(3) RIA=AI(1) RIB=AI(2) RIC=AI(3) AA=CONV3/RIA BB=CONV3/RIB CC=CONV3/RIC DO 8 K=1,3 RC(K)=CONV2/AI(K) 8 CONTINUE WRITE(IW,208)(RC(I),I=1,3),AA,BB,CC 208 FORMAT(//1x,78(1h-)/ * ' ROTATIONAL CONSTANTS A,B,C (cm-1 and MHz):'/ * 1x,78(1h-)/15X,3F20.8/15X,3F20.8) WRITE(IW,209) 209 FORMAT(//1x,78(1h-)/' INERTIAL DERIVATIVES: '/ * 1x,78(1h-)/9X,'J.aa',8X,'J.bb',8X, * 'J.cc',8X,'J.ab',8X,'J.ac',8X,'J.bc'/) WRITE(IW,210)(J,(A(I,J),I=1,6),J=1,NMOD) 210 FORMAT(1X,I2,'. ',6F12.7) C C C.............CENTRIFUGAL DISTORTION CONSTANTS.............. C C array TAU(i,j) is dimensioned for all 21 possible distinct TAU constants, C but at the moment only the following are defined: C TAU(i,1)=tau.aaaa TAU(i,2)=tau.bbbb TAU(i,3)=tau.cccc C TAU(i,4)=tau.aabb TAU(i,5)=tau.aacc TAU(i,6)=tau.bbcc C TAU(i,7)=tau.abab TAU(i,8)=tau.acac TAU(i,9)=tau.bcbc C TAU(i,10)=tau.bbbc TAU(i,11)=tau.abbb TAU(i,12)=tau.aaac C TAU(i,13)=tau.aaab TAU(i,14)=tau.accc TAU(i,15)=tau.bccc C where the first index gives the contribution to a given constant from C the i'th normal mode C RIA4=RIA**4 RIB4=RIB**4 RIC4=RIC**4 RIAB2=(RIA*RIB)**2 RIAC2=(RIA*RIC)**2 RIBC2=(RIB*RIC)**2 RIBBBC=RIB**3*RIC RIABBB=RIA*RIB**3 RIAAAC=RIA**3*RIC DO 12 I=1,NMOD TAU(I,1)=CONV4*(-0.5D0)*A(1,I)*A(1,I)/(RLAMBD(I,1)*RIA4) TAAAA=TAAAA+TAU(I,1) TAU(I,2)=CONV4*(-0.5D0)*A(2,I)*A(2,I)/(RLAMBD(I,1)*RIB4) TBBBB=TBBBB+TAU(I,2) TAU(I,3)=CONV4*(-0.5D0)*A(3,I)*A(3,I)/(RLAMBD(I,1)*RIC4) TCCCC=TCCCC+TAU(I,3) TAU(I,4)=CONV4*(-0.5D0)*A(1,I)*A(2,I)/(RLAMBD(I,1)*RIAB2) TAABB=TAABB+TAU(I,4) TAU(I,5)=CONV4*(-0.5D0)*A(1,I)*A(3,I)/(RLAMBD(I,1)*RIAC2) TAACC=TAACC+TAU(I,5) TAU(I,6)=CONV4*(-0.5D0)*A(2,I)*A(3,I)/(RLAMBD(I,1)*RIBC2) TBBCC=TBBCC+TAU(I,6) TAU(I,7)=CONV4*(-0.5D0)*A(4,I)*A(4,I)/(RLAMBD(I,1)*RIAB2) TABAB=TABAB+TAU(I,7) TAU(I,8)=CONV4*(-0.5D0)*A(5,I)*A(5,I)/(RLAMBD(I,1)*RIAC2) TACAC=TACAC+TAU(I,8) TAU(I,9)=CONV4*(-0.5D0)*A(6,I)*A(6,I)/(RLAMBD(I,1)*RIBC2) TBCBC=TBCBC+TAU(I,9) TAU(I,10)=CONV4*(-0.5D0)*A(2,I)*A(6,I)/(RLAMBD(I,1)*RIBBBC) TBBBC=TBBBC+TAU(I,10) TAU(I,11)=CONV4*(-0.5D0)*A(4,I)*A(2,I)/(RLAMBD(I,1)*RIABBB) TABBB=TABBB+TAU(I,11) TAU(I,12)=CONV4*(-0.5D0)*A(5,I)*A(1,I)/(RLAMBD(I,1)*RIAAAC) TAAAC=TAAAC+TAU(I,12) TAU(I,13)=CONV4*(-0.5D0)*A(1,I)*A(4,I)/(RLAMBD(I,1)*RIAAAC) TAAAB=TAAAB+TAU(I,13) TAU(I,14)=CONV4*(-0.5D0)*A(5,I)*A(3,I)/(RLAMBD(I,1)*RIAAAC) TACCC=TACCC+TAU(I,14) TAU(I,15)=CONV4*(-0.5D0)*A(6,I)*A(3,I)/(RLAMBD(I,1)*RIAAAC) TBCCC=TBCCC+TAU(I,15) 12 CONTINUE WRITE(IW,211) TAAAA*1000.D0,TBBBB*1000.D0,TCCCC*1000.D0, * TAABB*1000.D0,TAACC*1000.D0,TBBCC*1000.D0,TABAB*1000.D0, * TACAC*1000.D0,TBCBC*1000.D0, * TAAAB*1000.D0,TBBBC*1000.D0,TAAAC*1000.D0, * TABBB*1000.D0,TBCCC*1000.D0,TACCC*1000.D0 211 FORMAT(//1x,78(1h-)/' CENTRIFUGAL DISTORTION CONSTANTS /kHz:'/ * 1x,78(1h-)/ * ' T.aaaa, T.bbbb, T.cccc: ',7x,3F15.7/ * ' T.aabb, T.aacc, T.bbcc: ',7x,3F15.7/ * ' T.abab, T.acac, T.bcbc: ',7x,3F15.7/ * ' T.aaab, T.bbbc, T.aaac: ',7x,3F15.7/ * ' T.abbb, T.bccc, T.accc: ',7x,3F15.7) C C Watsonian constants: some elements of TAU are redefined to: C TAU(i,4)=tau'.aabb TAU(i,5)=tau'.aacc TAU(i,6)=tau'.bbcc C TAU(i,7)=tau.1 TAU(i,8)=tau.2 C such that C tau'.xxyy= tau.xxyy + 2 tau.xyxy, etc. C tau.1= tau'.yyzz + tau'.xxzz + tau'.xxyy C tau.2= B.x tau'.yyzz + B.y tau'.xxzz + B.z tau'.xxyy C [tau'.xxxx= tau.xxxx etc. ] C then C D.J=-(tau'.xxxx+tau'.yyyy)/8 C D.JK= (3 tau'.xxxx + 3 tau'.yyyy - 2 tau.1)/8 C D.K=-(tau'.xxxx + tau'.yyyy + tau'.zzzz - tau.1)/4 C delta.J=-(tau'.xxxx - tau'.yyyy)/16 C delta.K={ (B.x-B.z) tau'.xxxx + (B.y-B.z) tau'.yyyy - (B.x+B.y) tau.1 C + 2 tau.2 (B.x+B.y+B.z) } / { 8 (B.x-B.y) } C C These constants are calculated for two possible representations of the C Hamiltonian: C C prolate I.r: x <-> b oblate III.l: x <-> b C y <-> c y <-> a C z <-> a z <-> c C S=AA+BB+CC TAABBP=TAABB+2.D0*TABAB TAACCP=TAACC+2.D0*TACAC TBBCCP=TBBCC+2.D0*TBCBC TAU1=TAABBP+TAACCP+TBBCCP TAU2=(AA*TBBCCP+BB*TAACCP+CC*TAABBP)/S DO 14 I=1,NMOD TAU(I,4)=TAU(I,4)+2.D0*TAU(I,7) TAU(I,5)=TAU(I,5)+2.D0*TAU(I,8) TAU(I,6)=TAU(I,6)+2.D0*TAU(I,9) TAU(I,7)=TAU(I,4)+TAU(I,5)+TAU(I,6) TAU(I,8)=(AA*TAU(I,6)+BB*TAU(I,5)+CC*TAU(I,4))/S 14 CONTINUE WRITE(IW,212)TAABBP*1000.d0,TAACCP*1000.d0,TBBCCP*1000.d0, * TAU1*1000.d0,TAU2*1000.d0 212 FORMAT(/' T''.aabb, T''.aacc, T''.bbcc: ',5x,3F15.8/ * ' Tau.1, Tau.2 : ',5x,2F15.8/) C PDJ=-0.125D0*(TBBBB+TCCCC) PDJK=0.125D0*(3.D0*TBBBB+3.D0*TCCCC-2.D0*TAU1) PDK=-0.25D0*(TAAAA+TBBBB+TCCCC-TAU1) PDELJ=-0.0625D0*(TBBBB-TCCCC) IF(ABS((BB-CC)/BB).LT.1.D-5)THEN PDELK=0.0D0 GOTO 216 ENDIF PDELK=((BB-AA)*TBBBB+(CC-AA)*TCCCC-(BB+CC)*TAU1+2.D0*S*TAU2)/ * (8.D0*(BB-CC)) C ODJ=-0.125D0*(TBBBB+TAAAA) ODJK=0.125D0*(3.D0*TBBBB+3.D0*TAAAA-2.D0*TAU1) ODK=-0.25D0*(TAAAA+TBBBB+TCCCC-TAU1) ODELJ=-0.0625D0*(TBBBB-TAAAA) IF(ABS((BB-AA)/BB).LT.1.D-5)THEN ODELK=0.0D0 GOTO 216 ENDIF ODELK=((BB-CC)*TBBBB+(AA-CC)*TAAAA-(BB+AA)*TAU1+2.D0*S*TAU2)/ * (8.D0*(BB-AA)) C 216 ASYMK=(2.D0*BB-AA-CC)/(AA-CC) DO 16 I=1,NMOD IF(ASYMK.LE.0.D0)THEN DC(I,1)=-0.125D0*(TAU(I,2)+TAU(I,3)) DC(I,2)=0.125D0*(3.D0*TAU(I,2)+3.D0*TAU(I,3)-2.D0*TAU(I,7)) DC(I,3)=-0.25D0*(TAU(I,1)+TAU(I,2)+TAU(I,3)-TAU(I,7)) DC(I,4)=-0.0625D0*(TAU(I,2)-TAU(I,3)) IF(ABS(PDELJ).LT.1D-10)DC(I,4)=0.D0 IF(ABS((BB-CC)/BB).LT.1.D-10)THEN DC(I,5)=0.0D0 GOTO 16 ENDIF DC(I,5)=((BB-AA)*TAU(I,2)+(CC-AA)*TAU(I,3)-(BB+CC)*TAU(I,7)+ * 2.D0*S*TAU(I,8))/(8.D0*(BB-CC)) ELSE DC(I,1)=-0.125D0*(TAU(I,2)+TAU(I,1)) DC(I,2)=0.125D0*(3.D0*TAU(I,2)+3.D0*TAU(I,1)-2.D0*TAU(I,7)) DC(I,3)=-0.25D0*(TAU(I,1)+TAU(I,2)+TAU(I,3)-TAU(I,7)) DC(I,4)=-0.0625D0*(TAU(I,2)-TAU(I,1)) IF(ABS(ODELJ).LT.1D-10)DC(I,4)=0.D0 IF(ABS((BB-AA)/BB).LT.1.D-10)THEN DC(I,5)=0.0D0 GOTO 16 ENDIF DC(I,5)=((BB-CC)*TAU(I,2)+(AA-CC)*TAU(I,1)-(BB+AA)*TAU(I,7)+ * 2.D0*S*TAU(I,8))/(8.D0*(BB-AA)) ENDIF 16 CONTINUE C C...Optional input of experimental cd. constants C write(*,1001) 1001 format(1x/' QUARTIC CENTRIFUGAL DISTORTION CONSTANTS /kHz:') read(ir,'(a)',ERR=50,END=50)header if(header.ne.'CENTRIFUGA')THEN BACKSPACE(ir) goto 50 endif READ(IR,*,ERR=50,END=50)IASYM,(CDEXP(I),I=1,5) CDPRO(1)=PDJ CDPRO(2)=PDJK CDPRO(3)=PDK CDPRO(4)=PDELJ CDPRO(5)=PDELK CDOBL(1)=ODJ CDOBL(2)=ODJK CDOBL(3)=ODK CDOBL(4)=ODELJ CDOBL(5)=ODELK DO 80 I=1,5 ERRCD(I)=0.D0 80 CONTINUE C NCD=0 CHISQ=0.D0 DO 81 I=1,5 IF(CDEXP(I).NE.0.D0)THEN NCD=NCD+1 IF(IASYM.LE.0)THEN ERRCD(I)=(CDEXP(I)-CDPRO(I))/CDEXP(I) ELSE ERRCD(I)=(CDEXP(I)-CDOBL(I))/CDEXP(I) ENDIF CHISQ=CHISQ+ERRCD(I)**2 ENDIF 81 CONTINUE C sigma=dsqrt(chisq)/(NCD-1) C C...output of calculated and experimental cd. constants C IF(IASYM.LT.0)THEN WRITE(IW,51)'Prolate, I.r:',(CDEXP(I)*1000.d0,I=1,5), * (CDPRO(I)*1000.d0,I=1,5),(ERRCD(I),I=1,5),CHISQ,SIGMA WRITE(*,51)'Prolate, I.r:',(CDEXP(I)*1000.d0,I=1,5), * (CDPRO(I)*1000.d0,I=1,5),(ERRCD(I),I=1,5),CHISQ,SIGMA ELSE WRITE(IW,51)'Oblate, III.l:',(CDEXP(I)*1000.d0,I=1,5), * (CDPRO(I)*1000.d0,I=1,5),(ERRCD(I),I=1,5),CHISQ,SIGMA WRITE(*,51)'Oblate, III.l:',(CDEXP(I)*1000.d0,I=1,5), * (CDPRO(I)*1000.d0,I=1,5),(ERRCD(I),I=1,5),CHISQ,SIGMA ENDIF 51 FORMAT(1X,A/13X,'D.J ',9X,'D.JK',9X,'D.K ',9X,'d.J',9X,'d.K'/ * ' obs.',5F13.7/ * ' calc.',5F13.7/ * ' (o-c)/o',5F13.7/' chisq =',F13.8,' sigma =',f13.8) GOTO 53 C C...output of calculated cd. constants only C 50 WRITE(IW,214)PDJ*1000.d0,PDJK*1000.d0,PDK*1000.d0,PDELJ*1000.d0, * PDELK*1000.d0,ODJ*1000.d0,ODJK*1000.d0,ODK*1000.d0, * ODELJ*1000.d0,ODELK*1000.d0 WRITE(*,214)PDJ*1000.d0,PDJK*1000.d0,PDK*1000.d0,PDELJ*1000.d0, * PDELK*1000.d0,ODJ*1000.d0,ODJK*1000.d0,ODK*1000.d0, * ODELJ*1000.d0,ODELK*1000.d0 214 FORMAT(' ---------+'/ * ' Prolate | D.J, D.JK, DK : ',3F15.7/ * ' A, I.r | delta.J, delta.K: ',2F15.7/ * ' ---------+'/ * ' Oblate | D.J, D.JK, DK : ',3F15.7/ * ' A, III.l | delta.J, delta.K: ',2F15.7/ * ' ---------+' ) C C...S-reduction constants (only written to disk) c if(bb-cc.ne.0.d0)then psigma=(2.d0*aa-bb-cc)/(bb-cc) spd2=-pdelk/(4.d0*psigma) else spd2=0.d0 endif spd1=-pdelj spdj=pdj+2.d0*spd2 spdjk=pdjk-12.d0*spd2 spdk=pdk+10.d0*spd2 c if(bb-aa.ne.0.0d0)then osigma=(2.d0*cc-aa-bb)/(bb-aa) sod2=-odelk/(4.d0*osigma) else sod2=0.d0 endif sod1=-odelj sodj=odj+2.d0*sod2 sodjk=odjk-12.d0*sod2 sodk=odk+10.d0*sod2 WRITE(IW,1214)sPDJ*1000.d0,sPDJK*1000.d0,sPDK*1000.d0, * sPD1*1000.d0,sPD2*1000.d0,sODJ*1000.d0,sODJK*1000.d0, * sODK*1000.d0,sOD1*1000.d0,sOD2*1000.d0 1214 FORMAT(' ---------+'/ * ' Prolate | D.J, D.JK, DK : ',3F15.7/ * ' S, I.r | d1, d2 : ',2F15.7/ * ' ---------+'/ * ' Oblate | D.J, D.JK, DK : ',3F15.7/ * ' S, III.l | d1, d2 : ',2F15.7/ * ' ---------+' ) C C...output of vibrational mode contributions to cd. constants C 53 IF(ASYMK.LE.0.D0)THEN WRITE(IW,18)'(for prolate, I.r representation, constants)' ELSE WRITE(IW,18)'(for oblate, III.l representation, constants)' ENDIF 18 FORMAT(//1x,78(1h-)/ * ' VIBRATIONAL MODE CONTRIBUTIONS TO DISTORTION', * ' CONSTANTS (divided by value'/' of the constant):'/1X,A/ * 1x,78(1h-) /' mode',T16,'DJ',T30,'DJK',T44,'DK',T58,'dJ', * T72,'dK'/) IF(ASYMK.LE.0.0D0)THEN DJ=PDJ DJK=PDJK DK=PDK DELJ=PDELJ DELK=PDELK ELSE DJ=ODJ DJK=ODJK DK=ODK DELJ=ODELJ DELK=ODELK ENDIF IF(ABS(DELJ).LT.1.D-8)DELJ=1.D10 IF(ABS(DELK).LT.1.D-8)DELK=1.D10 DO 17 I=1,NMOD WRITE(IW,19)I,DC(I,1)/DJ,DC(I,2)/DJK,DC(I,3)/DK,DC(I,4)/DELJ, * DC(I,5)/DELK 17 CONTINUE 19 FORMAT(I3,'.',5X,5F14.8) C C C............CORIOLIS COUPLING CONSTANTS...................... C C zeta.a = l' M.a l (note that of the column eigenvectors in EIVR, C the first 5 or 6 are the rigid motions) C IA=1 IB=2 IC=3 NFMOD=N-NMOD+1 DO 30 I=1,3 DO 31 JA=NFMOD,N DO 31 JB=NFMOD,N ZSUM=0.d0 DO 32 JAT=1,NAT J1=3*(JAT-1)+IB J2=3*(JAT-1)+IC ZSUM=ZSUM-EIVR(J1,JA)*EIVR(J2,JB) 1 +EIVR(J2,JA)*EIVR(J1,JB) 32 CONTINUE ZETA(JA-NFMOD+1,JB-NFMOD+1,IA)=ZSUM 31 CONTINUE IF(IPRFUL.NE.0)THEN WRITE(IW,33)CHAR(ICHAR('a')+IA-1) 33 FORMAT(//1x,78(1h-)/' ZETA ',A,':'/1x,78(1h-)) c icors=1 icorf=nmod if(icorf.gt.14)icorf=14 960 WRITE(IW,133)(RLAMBD(JA,2),JA=icors,icorf) WRITE(IW,233)(JA,JA=icors,icorf) 233 FORMAT(12X,20(I7,1X)) 133 FORMAT(13X,20(F8.1 )) c DO 34 JA=1,NMOD do 962 jb=icors,icorf if(jb.eq.ja)then zeta(ja,jb,ia)=1000.D0 endif 962 continue c WRITE(IW,35)JA,RLAMBD(JA,2),(ZETA(JA,JB,IA),JB=icors,icorf) c do 963 jb=icors,icorf if(jb.eq.ja)then zeta(ja,jb,ia)=0.D0 endif 963 continue 34 CONTINUE 35 FORMAT(1X,I2,'.',F7.1,2x,20(F8.4)) if(icorf.lt.nmod)then icors=icors+14 icorf=icorf+14 if(icorf.gt.nmod)icorf=nmod write(iw,'(1x)') goto 960 endif ENDIF IA=IA+1 IB=IB+1 IC=IC+1 IF(IB.GT.3)IB=1 IF(IC.GT.3)IC=1 30 CONTINUE C C C............HARMONIC AND CORIOLIS CONTRIBUTIONS TO SEXTICS................ C C...evaluation of the cubic Coriolis terms C(k,alpha,alpha,alpha) C eq.(124),(103) [WATSON] - CAAAK is dimensionless, all terms are C converted first to SI C CAAA=0.D0 CBBB=0.D0 CCCC=0.D0 DO 300 I=1,3 DO 301 K=1,NMOD CAAAK(K,I)=0.D0 GAMMAK=DSQRT(RLAMBD(K,1)*6.0221366D+28)/1.05457266D-34 DO 302 L=1,NMOD GAMMAL=DSQRT(RLAMBD(L,1)*6.0221366D+28)/1.05457266D-34 CAAAK(K,I)=CAAAK(K,I)+(1.D0/GAMMAL**2+2.D0/GAMMAK**2) * *(A(I,L)*4.0749726D-24) *ZETA(L,K,I) 302 CONTINUE CAAAK(K,I)=CAAAK(K,I)*1.D0/( 4.D0*DSQRT(GAMMAK)* * (AI(I)*1.6605402D-47)**3) 301 CONTINUE 300 CONTINUE C C...Coriolis contribution to PHI(alpha,alpha,alpha), second term, first C equation in Table 7 [WATSON], first in units of OMEGA (cm-1), then MHz C phiac=0.D0 phibc=0.D0 phicc=0.D0 DO 303 K=1,NMOD PHIAC=PHIAC-2.D0*RLAMBD(K,2)*CAAAK(K,1)**2 PHIBC=PHIBC-2.D0*RLAMBD(K,2)*CAAAK(K,2)**2 PHICC=PHICC-2.D0*RLAMBD(K,2)*CAAAK(K,3)**2 303 CONTINUE PHIAC=PHIAC*29979.2458D0 PHIBC=PHIBC*29979.2458D0 PHICC=PHICC*29979.2458D0 c c...Harmonic contribution to PHI(alpha,alpha,alpha), first term, first C equation in Table 7 [WATSON] C PHIAH=(3.D0/16.D0)*( * TAAAA**2/AA+TAAAB**2/BB+TAAAC**2/CC) PHIBH=(3.D0/16.D0)*( * TABBB**2/AA+TBBBB**2/BB+TBBBC**2/CC) PHICH=(3.D0/16.D0)*( * TACCC**2/AA+TBCCC**2/BB+TCCCC**2/CC) c c...Contribution to PHI(alpha,alpha,alpha), from elimination of non- C orthorhombic terms, fourth term, first equation in Table 7,p.66 [WATSON] C if(aa.ne.cc)then PHIAR=0.25D0*(TAAAB**2/(AA-BB)+TAAAC**2/(AA-CC) ) else PHIAR=1.d10 endif if(bb.ne.cc)then PHIBR=0.25D0*(TABBB**2/(BB-AA)+TBBBC**2/(BB-CC) ) else PHIBR=1.D10 endif if(cc.ne.bb)then PHICR=0.25D0*(TACCC**2/(CC-AA)+TBCCC**2/(CC-BB) ) else PHICR=1.d10 endif C WRITE(IW,304)PHIAH*1000.D0,PHIBH*1000.d0,PHICH*1000.D0, * PHIAR*1000.D0,PHIBR*1000.D0,PHICR*1000.D0, * PHIAC*1000.D0,PHIBC*1000.D0,PHICC*1000.D0 304 FORMAT(//1x,78(1h-)/ * ' HARMONIC, TRANSFORMATION AND CORIOLIS CONTRIBUTIONS', * ' TO SEXTICS (kHz):'/1x,78(1h-)/ * 26x,'Phi.aaa',10x,'Phi.bbb',10x,'Phi.ccc'// * ' Harmonic:',3F17.11/ * ' Transformation:',3F17.11/ * ' Coriolis:',3F17.11) C C...When close to symmetric limit (h1 and h2 are small) then from Table 6, p.41 C [WATSON], or explicitly for a C.3v molecule eq.36 [ALWATS] C C H.J = (1/2) (PHI.xxx+PHI.yyy) C h.3 = (1/4) (PHI.xxx-PHI.yyy) C HJOBL=0.5d0*( (phiah+phiar+phiac)+(phibh+phibr+phibc) ) HJPRO=0.5d0*( (phibh+phibr+phibc)+(phich+phicr+phicc) ) H3OBL=0.25d0*( (phiah+phiar+phiac)-(phibh+phibr+phibc) ) H3PRO=0.25d0*( (phibh+phibr+phibc)-(phich+phicr+phicc) ) write(iw,305)(2*bb-aa-cc)/(aa-cc), * HJPRO*1.D9,HJOBL*1.D9,H3PRO*1.D9,H3OBL*1.D9 305 format(1x/' (kappa=',f8.5,')', * 26x,'prolate oblate'// * ' HJ /mHz (C3v symm, w/o anh. term)',2F17.8/ * ' h3 /mHz (C3v symm, w/o anh. term)',2F17.8) c C C...HARMONIC CONTRIBUTIONS TO MOMENTS OF INERTIA AND TO ROTATIONAL CONSTANTS c are defined by: c I* = I.v - Sum.s [(v.s + d.s/2)Eps.s] c Gordy&Cook 3rd ed., p.708. eq. 13.195 + 13.196 c for alpha.harm see: c Mills, in Molecular Spectroscopy: Modern Research, Academic Press, c London (1972): p.121, eq.13 c C NOTE: (d.Eps)harm give the *'rred moments of inertia, further addition of C (d.Eps)anh is necessary for I.e's. At the rotational constant C level alpha.harm is not a very useful pointer to the overall C alpha since it is generally smaller and can be of opposite C sign than alpha.anh. However inertial defect and associated planar c moment are evaluated reliably for near-planar molecules since in c such case anharmonic contributions cancel. C DO 41 I=1,6 TEMP(I)=0.D0 41 CONTINUE DO 36 IA=1,3 DO 37 JA=1,NMOD ASUM=0.D0 DO 38 I=1,3 IF(IA.EQ.I)INDEX=I IF(IA.NE.I)INDEX=1+IA+I ASUM=ASUM+A(INDEX,JA)**2/AI(I) 38 CONTINUE ASUM=0.75D0*ASUM ZSUM=0.D0 DO 39 JB=1,NMOD IF(JA.EQ.JB)GOTO 39 IF(RLAMBD(JA,1).EQ.RLAMBD(JB,1))GOTO 39 ZSUM=ZSUM+ZETA(JA,JB,IA)**2 * *(3.D0*RLAMBD(JA,1)+RLAMBD(JB,1)) * /(RLAMBD(JA,1)-RLAMBD(JB,1)) 39 CONTINUE DEPSIL(JA,IA)=-2.D0*RK*(ASUM+ZSUM)/RLAMBD(JA,2) ALPHA(JA,IA)=DEPSIL(JA,IA)*CONV3/AI(IA)**2 TEMP(IA)=TEMP(IA)+DEPSIL(JA,IA) TEMP(IA+3)=TEMP(IA+3)+ALPHA(JA,IA) 37 CONTINUE 36 CONTINUE C WRITE(IW,40) 40 FORMAT(//1x,78(1h-)/ * ' VIBRATIONAL-ROTATIONAL CONTRIBUTIONS (d.Eps)harm ', * 'and (alpha)harm'/ * ' (in u Angstr^2 AND MHz resp.) defined by:'/ * ' I.v = I.e + Sum_s [( v.s + d.s/2)Eps.s ] and'/ * ' B.v = B.e - Sum_s [( v.s + d.s/2)alpha.s ]'/, * 1x,78(1h-)/ * ' MODE',7X,'eps.a',5X,'eps.b',5X,'eps.c',7x,'dv'8x, * 'alpha.A',2X,'alpha.B',2X,'alpha.C cm-1 dv+Delta_0'/) dzero=0.5d0*(temp(3)-temp(2)-temp(1)) DO 42 JA=1,NMOD dv= (depsil(ja,3)-depsil(ja,2)-depsil(ja,1)) WRITE(IW,43)JA,(DEPSIL(JA,I),I=1,3),dv,(ALPHA(JA,I),I=1,3), * rlambd(ja,2),dv+dzero 42 CONTINUE 43 FORMAT(1X,I2,'. ',3F10.4,F12.4,3X,3F9.2,f12.2,f12.5) dv= (temp(3)-temp(2)-temp(1)) WRITE(IW,44)(TEMP(I),I=1,3),dv,(TEMP(I),I=4,6) 44 FORMAT(1X/' Total:',3F10.4,F12.4,3x,3F9.2/) DO 45 I=1,3 SGN='+' ZSUM=-0.5D0*TEMP(I) IF(ZSUM.LT.0.D0)SGN='-' WRITE(IW,46)CHAR(ICHAR('a')+I-1),CHAR(ICHAR('a')+I-1),SGN, * DABS(ZSUM) 45 CONTINUE 46 FORMAT(1X,'(I.',A,')^* = (I.',A,')_0 ',A,F14.6, * ' u Angstr^2') c write(IW,47)dzero 47 format(1x/' Delta_0 =',16x,f14.6,' u Angstr^2') C RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________ C C GEOMAS - geometry and masses C VIBFR - vibrational secular equations C TWIST - convert general torsional coordinate to Cartesian displacements C STWIST - convert simple torsional coordinate to Cartesian dispalcements C DIAG - matrix diagonalisation C ROTCS - rotational transformation matrix C DCS - unit vectors along bond defining an internal coordinate C CROSS - vector product of unit vectors C DOT - sclar product of unit vectors (function) C MMULT - multiply two matrices C VRCON - various vibration-rotation constants C_____________________________________________________________________________