C------------------------------------------------------------------------------- C C CORSCL - scaling of molecular geometry by a global scaling factor C C------------------------------------------------------------------------------- C C The main purpose of this program is to enhance the accuracy of prediction C of rotational spectra for isotopic species of a molecule before its C precise geometry has been determined. C C The program allows scaling of ab initio, ar any other, trial geometry C for the molecule to best reproduce experimental rotational constants. C C C ver. 8a.IX.1998 ----- Zbigniew KISIEL ----- C ----- Marek JAWORSKI ----- 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 C Two input files are necessary: C C MOLNAM.COR - user defined file containing the Cartesians and masses of C atoms in the molecule in the .COR standard used by PMIFST. C C CORSCL.INP - single line ASCII file containing the values of rotational C constants A,B,C (in MHz) for the declared isotopic species. C C C The structure scaling factor can then be chosen manually for best C reproduction of the experimental moments of inertia (defined C the declared rotational constants) and the scaled Cartesians are C written to file CORSCL.OUT for subsequent use in PMIFST C C C Once the value of the scaling factor is established the output of scaled C data is made to file CORSCL.OUT C C------------------------------------------------------------------------------- C If program crashes on input then use the coordinates as printed out C into PMIFST.RES C C If in doubt use the yardstick below for aligning NATOM,X,Y,Z,WT. C The NATOM label is an arbitrary four digit descriptive string. C C...////////////////................///////////////.................... C C C------------------------------------------------------------------------------- C implicit real*8 (a-h,o-z) c parameter (maxat=100,conv=505379.07d0) real*8 cordin(maxat,4),refpmi(3,3) character*79 title character*30 fname1 character*4 latom(maxat) COMMON /MOMIA/CARTDF(maxat,4),pmi(3,3,2) WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | CORSCL - Scaling of cartesians in a .COR structural ', * T79,'|'/ * ' | file by a global scaling factor', * T79,'|'/ * ' |',76(1H_),'|'/' version 8a.IX.1998',T64,'Zbigniew KISIEL'/ * T64,'Marek JAWORSKI'//) 105 write(*,100)' Name of .COR data file: ' 100 format(1x,a,$) read(*,101)fname1 open(1,file=fname1,status='old',err=105) c open(3,file='corscl.inp',status='old',err=21) read(3,*)(refpmi(i,2),i=1,3) do 22 i=1,3 refpmi(i,1)=conv/refpmi(i,2) 22 continue close(3) c c c...read the .COR file c read(1,101)title 101 format(a) read(1,102)natoms 102 format(i4) do 20 i=1,natoms read(1,104)latom(i),(cordin(i,j),j=1,4) 20 continue 104 FORMAT(A4,3F16.10,f20.10) close(1) c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c...scaling c 11 write(*,'(1x/)') 12 write(*,100)' scaling factor (ENTER to end) = ' read(*,'(1f20.10)',err=12)scale if(scale.eq.0.0)goto 10 sclast=scale c do 14 i=1,natoms cartdf(i,4)=cordin(i,4) do 14 j=1,3 cartdf(i,j)=scale*cordin(i,j) 14 continue c call momina(natoms,1) write(*,23)'Rotational constants:', * (refpmi(i,2),i=1,3),(pmi(i,2,1),i=1,3), * (refpmi(i,2)-pmi(i,2,1),i=1,3) write(*,23)'Moments of inertia:', * (refpmi(i,1),i=1,3),(pmi(i,1,1),i=1,3), * (refpmi(i,1)-pmi(i,1,1),i=1,3) 23 format(1x/1x,a/' Ref: ',3f15.6/ * ' Calc: ',3f15.6/ * ' Ref-Calc: ',3f15.6) sigma=0.0d0 do 24 n=1,3 sigma=sigma+(refpmi(n,1)-pmi(n,1,1))**2 24 continue write(*,25)sigma 25 format(1x/' SUM I(ref-calc)**2 = ',f20.6) c goto 11 c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...output to file and exit c 10 open(2,file='corscl.out',status='unknown') c write(2,101)title write(2,102)natoms do 15 n=1,natoms write(2,103)latom(n),(cartdf(n,j),j=1,4) 15 continue 103 format(a4,3(2x,f13.10),f14.8) write(2,'('' |''/'' |__ Scaled with scaling factor = '', * f10.7/)')sclast c write(2,23)'Rotational constants:', * (refpmi(i,2),i=1,3),(pmi(i,2,1),i=1,3), * (refpmi(i,2)-pmi(i,2,1),i=1,3) write(2,23)'Moments of inertia:', * (refpmi(i,1),i=1,3),(pmi(i,1,1),i=1,3), * (refpmi(i,1)-pmi(i,1,1),i=1,3) write(2,25)sigma c close(2) write(*,'(1x//)') stop c 21 write(*,'(// * '' **** ERROR: Nonexistant or problematic CORSCL.INP file''//)') c stop end C_____________________________________________________________________________ C SUBROUTINE MOMINA(natms,icalc) C C Principal moments of inertia of a molecule from its Cartesians and c masses supplied via CARTDF in the COMMON block /MOMIA/ C C natms = number of atoms c icalc = 1 or 2 defines the output data set (last index of PMI) c PMI = first index defines a,b,c c = second index defines moment of inertia, rotational constant, c or planar moment c = third index defines number of data set C parameter (maxat=100) IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MOMIA/CARTDF(maxat,4),pmi(3,3,2) REAL*8 CORD(maxat,3),A(6),EIGVEC(3,3) CON=0.1745329252199D-1 C C...Transformation to centre of mass coordinates C 6000 WW=0.0D0 DO 991 I=1,NATMS WW=WW+CARTDF(I,4) cord(i,1)=cartdf(i,1) cord(i,2)=cartdf(i,2) cord(i,3)=cartdf(i,3) 991 CONTINUE FMX=0.0D0 FMY=0.0D0 FMZ=0.0D0 DO 992 N=1,NATMS FMX=FMX+CARTDF(N,4)*CORD(N,1) FMY=FMY+CARTDF(N,4)*CORD(N,2) FMZ=FMZ+CARTDF(N,4)*CORD(N,3) 992 CONTINUE DO 993 N=1,NATMS CORD(N,1)=CORD(N,1)-FMX/WW CORD(N,2)=CORD(N,2)-FMY/WW CORD(N,3)=CORD(N,3)-FMZ/WW 993 CONTINUE C C Principal moments of inertia, rotational constants and planar moments C DO 994 N=1,6 A(N)=0.0D0 994 CONTINUE DO 995 N=1,NATMS A(1)=A(1)+CARTDF(N,4)*(CORD(N,2)**2+CORD(N,3)**2) A(2)=A(2)-CARTDF(N,4)*CORD(N,1)*CORD(N,2) A(3)=A(3)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,3)**2) A(4)=A(4)-CARTDF(N,4)*CORD(N,1)*CORD(N,3) A(5)=A(5)-CARTDF(N,4)*CORD(N,2)*CORD(N,3) A(6)=A(6)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,2)**2) 995 CONTINUE C CALL EIGEN(A,EIGVEC,3,0) C PMI(1,1,icalc)=A(6) PMI(2,1,icalc)=A(3) PMI(3,1,icalc)=A(1) PMI(1,2,icalc)=505379.07D0/PMI(1,1,icalc) PMI(2,2,icalc)=505379.07D0/PMI(2,1,icalc) PMI(3,2,icalc)=505379.07D0/PMI(3,1,icalc) PMI(1,3,icalc)=0.5d0* 1 (pmi(2,1,icalc)+pmi(3,1,icalc)-pmi(1,1,icalc)) PMI(2,3,icalc)=0.5d0* 1 (pmi(1,1,icalc)+pmi(3,1,icalc)-pmi(2,1,icalc)) PMI(3,3,icalc)=0.5d0* 1 (pmi(1,1,icalc)+pmi(2,1,icalc)-pmi(3,1,icalc)) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE EIGEN(A,R,N,MV) C C MATRIX DIAGONALIZATION ROUTINE FOR SYMMETRIC MATRICES FROM THE C 'SSP' LIBRARY C C - Jacobi type algorithm C - matrix to be diagonalized has to be in symmetric storage mode C (A and R are presently dimensioned for 10x10 matrices) C DIMENSION A(55),R(100) DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, 1 COSX2,SINCS,RANGE,DSQRT,DABS C C C GENERATE IDENTITY MATRIX C 5 RANGE=1.0D-12 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0d0 IF(I-J) 20,15,20 15 R(IJ)=1.0d0 20 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) C 25 ANORM=0.0d0 DO 35 I=1,N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM) 165,165,40 40 ANORM=1.414213562d0*DSQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR C IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 C C COMPUTE SIN AND COS C 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(DABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=0.5d0*(A(LL)-A(MM)) 68 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/DSQRT(2.0d0*(1.0d0+(DSQRT(1.0d0-Y*Y)))) SINX2=SINX*SINX 78 COSX=DSQRT(1.0d0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.0d0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLETION C C TEST FOR M = LAST COLUMN C 130 IF(M-N) 135,140,135 135 M=M+1 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C 160 IF(THR-ANRMX) 165,165,45 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE C RETURN END C_____________________________________________________________________________ C_____________________________________________________________________________