C------------------------------------------------------------------------------- C C CORSCL - scaling of molecular geometry 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 Provision for automatic generation of constants for several types of C isotopically substituted species. C C C ver. 17.V.2017 ----- 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 principal coordinates C and masses of atoms in the molecule in the .COR standard C 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 as follows: C C 1/ manually derived single scaling factor for best reproduction of C the experimental moments of inertia (as defined by the rotational C constants in CORSCL.INP) C 2/ automatic three axis scaling based on planar moments making calculated C moments of inertai identically equal to the reference ones from CORSCL.INP C C C The scaled a,b,c are written to file with the input name extended by '_scaled.out'. C C------------------------------------------------------------------------------- C C Modification history: C C ????: Creation, mid 90's? C 24.10.12: Major revision including option of exact scaling using planar moments C and predictions for single isotopic substitutions C 7.03.13: More isotopic substitution options C 9.08.14: Incremental improvements C 23.02.17: Debugging of crash and modified calc. for planar molecules c 17.05.17: More informative input error diagnostics C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C implicit real*8 (a-h,o-z) c parameter (maxat=100,conv=505379.009d0,maxiso=500) real*8 cordin(maxat,4),refpmi(3,3),unspmi(3,3) real*8 acon(maxiso),bcon(maxiso),ccon(maxiso), * pa(maxiso),pb(maxiso),pc(maxiso) character title*79,line*100,fname1*100,fname2*120,isonam*2 character errmes*80 character*4 latom(maxat) integer numat(maxiso),numata(maxiso) COMMON /MOMIA/CARTDF(maxat,4),pmi(3,3) c scale=1.0d0 deltai=0.d0 c c...header c WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | CORSCL - Scaling of Cartesians from a .COR structural ', * T79,'|'/ * ' | file of atomic coordinates', * T79,'|'/ * ' |',76(1H_),'|'/' version 17.V.2017',T64,'Zbigniew KISIEL'/ * T64,'Marek JAWORSKI'//) c write(*,63) 63 format(' NOTES: ', * 'The input Cartesians can be in arbitrary axes but'/ * 8x,'the output will be in principal coordinates.'// * 8X,'The first four characters of each atom declaration line'/ * 8x,'are reserved for an arbitrary alphanumeric descriptor'///) c 105 write(*,100)'Name of the data file (.COR standard of PMIFST): ' 100 format(1x,a,$) read(*,101)fname1 open(1,file=fname1,status='old',err=105) c c c...read the experimental A,B,C c open(3,file='corscl.inp',status='old',err=21) read(3,*,err=221)(refpmi(i,2),i=1,3) abc=refpmi(1,2)*refpmi(2,2)*refpmi(3,2) if(abc.le.0.d0)goto 521 c do 22 i=1,3 refpmi(i,1)=conv/refpmi(i,2) 22 continue refpmi(1,3)=0.5d0*(refpmi(2,1)+refpmi(3,1)-refpmi(1,1)) refpmi(2,3)=0.5d0*(refpmi(1,1)+refpmi(3,1)-refpmi(2,1)) refpmi(3,3)=0.5d0*(refpmi(1,1)+refpmi(2,1)-refpmi(3,1)) close(3) c c c...read the .COR file c errmes='Reading the title from the .COR file' read(1,101,err=321,end=321)title 101 format(a) errmes='Reading the no of atoms from the .COR file' read(1,*,err=321,end=321)natoms do 20 i=1,natoms errmes='Reading an atom definition line from the .COR file' read(1,'(a)',err=321,end=421)line latom(i)=line(1:4) line=line(5:) errmes='Reading from an atom definition line from the .COR file' read(line,*,err=321,end=421)(cordin(i,j),j=1,4) 20 continue 104 FORMAT(A4,1x,4F16.10) close(1) ifirst=1 fname2=fname1(1:len_trim(fname1))//'_scaled.out' c c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c MAIN LOOP c c...options c 11 write(*,'(1x/)') if(ifirst.eq.1)then write(*,100)' scaling factor = 1.' scale=1.d0 ifirst=0 else 12 write(*,1000) 1000 format( * ' Type one of: '//8x, * ' x.xx = global scaling factor x.xx to multiply the input'/ * 19x, 'coordinates'/8x, * ' 0 = just exit'/8x, * ' -1 = save input coordinates scaled with the last'/19x, * 'scaling factor'/8x, * ' -2 = scale input coordinates using planar moments and'/ * 19x, 'calculate isotopic species'//8x * '..... ',$) read(*,*,err=12)scale c c...exit to save options c if(nint(scale).eq.-1)goto 10 if(nint(scale).eq.-2)goto 50 if(scale.eq.0.d0)stop if(scale.lt.0.d0)goto 12 endif c c...new scaling c 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) <----- c if(scale.eq.1.0d0)then do 60 i=1,3 do 60 j=1,3 unspmi(i,j)=pmi(i,j) 60 continue do 65 i=1,natoms do 65 j=1,3 cordin(i,j)=cartdf(i,j) 65 continue endif c write(*,23)'Rotational constants:', * (unspmi(i,2),i=1,3), * (refpmi(i,2),i=1,3), * (pmi(i,2),i=1,3), * (refpmi(i,2)-pmi(i,2),i=1,3) write(*,23)'Moments of inertia:', * (unspmi(i,1),i=1,3), * (refpmi(i,1),i=1,3), * (pmi(i,1),i=1,3), * (refpmi(i,1)-pmi(i,1),i=1,3) write(*,23)'Planar moments:', * (unspmi(i,3),i=1,3), * (refpmi(i,3),i=1,3), * (pmi(i,3),i=1,3), * (refpmi(i,3)-pmi(i,3),i=1,3) c c...Scaling coefficients based on planar moments. c c If a planar moment is zero then the reference inertial defect is c used in the calculation of Ic c if(scale.eq.1.0d0)then if(dabs(pmi(1,3)).lt.1.d-10)then scalea=0.d0 else scalea=dsqrt(refpmi(1,3)/pmi(1,3)) endif c if(dabs(pmi(2,3)).le.1.d-10)then scaleb=0.d0 else scaleb=dsqrt(refpmi(2,3)/pmi(2,3)) endif c if(dabs(pmi(3,3)).le.1.d-10)then scalec=0.d0 deltai=(refpmi(3,1)-refpmi(1,1)-refpmi(2,1)) c * *(scalea+scaleb)/2.d0 else scalec=dsqrt(refpmi(3,3)/pmi(3,3)) endif c write(*,'(1x/'' possible axis scaling = ''/ * 15x,3f15.6)')scalea,scaleb,scalec endif c 23 format(1x/1x,a/' Unscaled: ',3F15.6/ * ' Ref: ',3f15.6/ * ' Scaled: ',3f15.6/ * ' Ref-Scaled: ',3f15.6) sigma=0.0d0 do 24 n=1,3 sigma=sigma+(refpmi(n,1)-pmi(n,1))**2 24 continue write(*,25)sigma 25 format(1x/' SUM I(ref-scld)**2 = ',f20.6) c goto 11 c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c c...SCALE=-1: output of globally scaled coordinates to the .OUT file and exit c 10 open(2,file=fname2(1:len_trim(fname2)),status='unknown') c write(2,101)title write(2,102)natoms 102 format(i5) 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:', * (unspmi(i,2),i=1,3), * (refpmi(i,2),i=1,3), * (pmi(i,2),i=1,3), * (refpmi(i,2)-pmi(i,2),i=1,3) write(2,23)'Moments of inertia:', * (unspmi(i,1),i=1,3), * (refpmi(i,1),i=1,3), * (pmi(i,1),i=1,3), * (refpmi(i,1)-pmi(i,1),i=1,3) write(2,23)'Planar moments:', * (unspmi(i,3),i=1,3), * (refpmi(i,3),i=1,3), * (pmi(i,3),i=1,3), * (refpmi(i,3)-pmi(i,3),i=1,3) write(2,25)sigma c close(2) write(*,'(1x//1x,2a)')'Output written to: ', * fname2(1:len_trim(fname2)) write(*,'(1x/)') stop c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c c...SCALE=-2: output of planar moment scaled coordinates to the .OUT file c 50 do 54 i=1,natoms cartdf(i,4)=cordin(i,4) cartdf(i,1)=scalea*cordin(i,1) cartdf(i,2)=scaleb*cordin(i,2) cartdf(i,3)=scalec*cordin(i,3) 54 continue c open(2,file=fname2(1:len_trim(fname2)),status='unknown') c write(2,101)title write(2,102)natoms do 55 n=1,natoms write(2,103)latom(n),(cartdf(n,j),j=1,4) 55 continue write(2,155)scalea,scaleb,scalec 155 format(' |'/' |__ Scaled with scaling factors from ', * 'planar moments:'// * 5x,'a-axis = ',f16.10/ * 5x,'b-axis = ',f16.10/ * 5x,'c-axis = ',f16.10) c call momina(natoms) <----- c if(scalec.eq.0.d0)then if planar pmi(3,3)=-0.5*deltai Pc c pmi(1,1)=pmi(2,3)+pmi(3,3) Ia=(Pb+Pc) pmi(1,2)=CONV/pmi(1,1) A pmi(2,1)=(pmi(1,3)+pmi(3,3)) Ib pmi(2,2)=CONV/pmi(2,1) B c pmi(3,1)=pmi(1,1)+pmi(2,1)+deltai Ic pmi(3,2)=CONV/pmi(3,1) C endif c write(2,23)'Rotational constants:', * (unspmi(i,2),i=1,3), * (refpmi(i,2),i=1,3), * (pmi(i,2),i=1,3), * (refpmi(i,2)-pmi(i,2),i=1,3) write(2,23)'Moments of inertia:', * (unspmi(i,1),i=1,3), * (refpmi(i,1),i=1,3), * (pmi(i,1),i=1,3), * (refpmi(i,1)-pmi(i,1),i=1,3) write(2,23)'Planar moments:', * (unspmi(i,3),i=1,3), * (refpmi(i,3),i=1,3), * (pmi(i,3),i=1,3), * (refpmi(i,3)-pmi(i,3),i=1,3) c sigma=0.d0 do 61 n=1,3 sigma=sigma+(refpmi(n,1)-pmi(n,1))**2 61 continue write(2,25)sigma c c c...save parent species constants c nisot=1 numat(nisot)=0 ACON(nisot)=pmi(1,2) BCON(nisot)=pmi(2,2) CCON(nisot)=pmi(3,2) PA(nisot)=pmi(1,3) PB(nisot)=pmi(2,3) PC(nisot)=pmi(3,3) c c...decide on isotopic substitution c 1002 write(*,1001) 1001 format(1x/' Type one of:'// * ' 0 or ENTER = output for parent isotopic species'/ * ' 2 = and all single 2H species'/ * ' 13 = and all single 13C species'/ * ' 15 = and all single 15N species'/ * ' 16 = and all single 16O species', * ' (for 18O parent)'/ * ' 17 = and all single 17O species'/ * ' 18 = and all single 18O species'/ * ' 37 = and all single 37Cl species'// * ' (a negative value will generate all'/ * ' double isotopic species)'// * ' ..... ',$) read(*,'(i5)',err=1002)isocon c isubst=1 if(isocon.lt.0)isubst=2 c if(iabs(isocon).eq.2)then rmpar=1.00782503207D0 mpar=1 rmiso=2.0141017778D0 isonam='H ' goto 1006 endif c if(iabs(isocon).eq.13)then rmpar=12.D0 mpar=12 rmiso=13.0033548378D0 isonam='C ' goto 1006 endif c if(iabs(isocon).eq.15)then rmpar=14.0030740052D0 mpar=14 rmiso=15.0001088984D0 isonam='N ' goto 1006 endif c if(iabs(isocon).eq.16)then rmpar=17.9991604D0 mpar=18 rmiso=15.9949146221D0 isonam='O ' goto 1006 endif c if(iabs(isocon).eq.17)then rmpar=15.9949146221D0 mpar=17 rmiso=16.99913150D0 isonam='O ' goto 1006 endif c if(iabs(isocon).eq.18)then rmpar=15.9949146221D0 mpar=16 rmiso=17.9991604D0 isonam='O ' goto 1006 endif c if(iabs(isocon).eq.37)then rmpar=34.96885268D0 mpar=35 rmiso=36.96590259D0 isonam='Cl' goto 1006 endif c if(isocon.eq.0)goto 1003 goto 1002 c c...optional constants for all singly substituted species for the selected c isotope c 1006 if(isubst.eq.1)then do 500 i=1,natoms if(nint(cartdf(i,4)).ne.dble(mpar))goto 500 cartdf(i,4)=rmiso nisot=nisot+1 call momina(natoms) <----- numat(nisot)=i c if(scalec.eq.0.d0)then if planar pmi(3,3)=-0.5*deltai Pc c pmi(1,1)=pmi(2,3)+pmi(3,3) Ia=(Pb+Pc) pmi(1,2)=CONV/pmi(1,1) A pmi(2,1)=(pmi(1,3)+pmi(3,3)) Ib pmi(2,2)=CONV/pmi(2,1) B c pmi(3,1)=pmi(1,1)+pmi(2,1)+deltai Ic pmi(3,2)=CONV/pmi(3,1) C endif c ACON(nisot)=pmi(1,2) BCON(nisot)=pmi(2,2) CCON(nisot)=pmi(3,2) PA(nisot)=pmi(1,3) PB(nisot)=pmi(2,3) PC(nisot)=pmi(3,3) cartdf(i,4)=rmpar 500 continue endif c c...optional constants for all doubly substituted species for the selected c isotope c if(isubst.eq.2)then do 501 i=1,natoms-1 if(nint(cartdf(i,4)).ne.dble(mpar))goto 501 do 510 j=i+1,natoms if(nint(cartdf(j,4)).ne.dble(mpar))goto 510 cartdf(i,4)=rmiso cartdf(j,4)=rmiso nisot=nisot+1 call momina(natoms) <----- numat(nisot)=i numata(nisot)=j c if(scalec.eq.0.d0)then if planar pmi(3,3)=-0.5*deltai Pc c pmi(1,1)=pmi(2,3)+pmi(3,3) Ia=(Pb+Pc) pmi(1,2)=CONV/pmi(1,1) A pmi(2,1)=(pmi(1,3)+pmi(3,3)) Ib pmi(2,2)=CONV/pmi(2,1) B c pmi(3,1)=pmi(1,1)+pmi(2,1)+deltai Ic pmi(3,2)=CONV/pmi(3,1) C endif c ACON(nisot)=pmi(1,2) BCON(nisot)=pmi(2,2) CCON(nisot)=pmi(3,2) PA(nisot)=pmi(1,3) PB(nisot)=pmi(2,3) PC(nisot)=pmi(3,3) cartdf(i,4)=rmpar cartdf(j,4)=rmpar 510 continue 501 continue endif c c c...output c 1003 write(2,532) 532 format(1x/80(1H-)) if(isocon.eq.0)then write(2,504)fname1(1:len_trim(fname1)) 504 format(' Scaled prediction for the parent ', * 'species using data'/' from file: ',a) else write(2,'(1x)') if(isocon.gt.1)then write(2,505)isocon,isonam,fname1(1:len_trim(fname1)) 505 format(' Scaled prediction for the parent and all single', * i3,a,' species using data'/' from file: ',a) endif c if(isocon.lt.0)then write(2,525)iabs(isocon),isonam,fname1(1:len_trim(fname1)) 525 format(' Scaled prediction for the parent and all double', * i3,a,' species using data'/' from file: ',a) endif endif c nisost=1 if(nisot.lt.10)then nisoen=nisot else nisoen=10 endif 511 if(nisot.eq.1)then write(2,502)(numat(i),i=1,nisot) else write(2,'(1x/'' Rotational constants and planar moments:'')') write(2,1502)iabs(isocon),isonam,(numat(i),i=nisost,nisoen) if(isubst.eq.2) * write(2,1502)iabs(isocon),isonam,(numata(i),i=nisost,nisoen) endif write(2,'(1x)') write(2,503)' A ',(acon(i),i=nisost,nisoen) write(2,503)' B ',(bcon(i),i=nisost,nisoen) write(2,503)' C ',(ccon(i),i=nisost,nisoen) write(2,'(1x)') write(2,503)' Pa ',(pa(i),i=nisost,nisoen) write(2,503)' Pb ',(pb(i),i=nisost,nisoen) write(2,503)' Pc ',(pc(i),i=nisost,nisoen) c if(iabs(isocon).gt.1)then write(2,'(1x/'' Differences from parent:'')') write(2,1502)iabs(isocon),isonam,(numat(i),i=nisost,nisoen) if(isubst.eq.2) * write(2,1502)iabs(isocon),isonam,(numata(i),i=nisost,nisoen) write(2,'(1x)') write(2,503)' A ',(acon(i)-acon(1),i=nisost,nisoen) write(2,503)' B ',(bcon(i)-bcon(1),i=nisost,nisoen) write(2,503)' C ',(ccon(i)-ccon(1),i=nisost,nisoen) write(2,'(1x)') write(2,503)' Pa ',(pa(i)-pa(1),i=nisost,nisoen) write(2,503)' Pb ',(pb(i)-pb(1),i=nisost,nisoen) write(2,503)' Pc ',(pc(i)-pc(1),i=nisost,nisoen) endif c if(nisoen.lt.nisot)then nisost=nisost+10 nisoen=nisoen+10 if(nisoen.gt.nisot)nisoen=nisot goto 511 endif c 502 format(1x/4x,50(i6,6x)) 1502 format(i2,a2,50(i6,6x)) 503 format(a,50f12.5) c if(nisot.gt.1)then write(*,530)nisot-1 write(2,531)nisot-1 530 format(1x/' Total number of predicted substituted species =',i5) 531 format(1x/' Total number of predicted substituted species =',i5/ * ' Atoms have been labeled according to their order in the ', * 'input file'//80(1H-)/) endif close(2) c write(*,'(1x//1x,2a)')'Output written to: ', * fname2(1:len_trim(fname2)) write(*,'(1x/)') stop c c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c 21 write(*,'(// * '' **** ERROR: Nonexistant or inaccessible CORSCL.INP file''//)') stop c 221 write(*,'(// * '' **** ERROR: Cannot read A,B,C from the CORSCL.INP file''//)') stop c 321 write(*,'(//'' **** ERROR: '',a//1x,a//)') * errmes(1:len_trim(errmes)),line(1:len_trim(line)) stop c 421 write(*,'(// * '' **** ERROR: Premature end of the coordinates file''/ * '' when reading atom number'',i5//)')i stop c 521 write(*,'(// *'' **** ERROR: Zero or negative A,B, or C in the CORSCL.INP file'' * //)') stop c end C C_____________________________________________________________________________ C SUBROUTINE MOMINA(natms) 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 PMI = first index defines a,b,c c = second index defines moment of inertia, rotational constant, c or planar moment C parameter (maxat=100) IMPLICIT REAL*8 (A-H,O-Z) C COMMON /MOMIA/CARTDF(maxat,4),pmi(3,3) REAL*8 CORD(maxat,3),A(6),EIGVEC(9),EIGINV(3,3), * cordx(maxat),cordy(maxat),cordz(maxat) 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)=A(6) PMI(2,1)=A(3) PMI(3,1)=A(1) PMI(1,2)=505379.009D0/PMI(1,1) PMI(2,2)=505379.009D0/PMI(2,1) PMI(3,2)=505379.009D0/PMI(3,1) PMI(1,3)=0.5d0* 1 (pmi(2,1)+pmi(3,1)-pmi(1,1)) PMI(2,3)=0.5d0* 1 (pmi(1,1)+pmi(3,1)-pmi(2,1)) PMI(3,3)=0.5d0* 1 (pmi(1,1)+pmi(2,1)-pmi(3,1)) C C Principal coordinates C EIGINV(1,1)=EIGVEC(5)*EIGVEC(3)-EIGVEC(6)*EIGVEC(2) EIGINV(1,2)=EIGVEC(6)*EIGVEC(1)-EIGVEC(4)*EIGVEC(3) EIGINV(1,3)=EIGVEC(4)*EIGVEC(2)-EIGVEC(5)*EIGVEC(1) EIGINV(2,1)=EIGVEC(2)*EIGVEC(9)-EIGVEC(8)*EIGVEC(3) EIGINV(2,2)=EIGVEC(7)*EIGVEC(3)-EIGVEC(9)*EIGVEC(1) EIGINV(2,3)=EIGVEC(8)*EIGVEC(1)-EIGVEC(7)*EIGVEC(2) EIGINV(3,1)=EIGVEC(8)*EIGVEC(6)-EIGVEC(5)*EIGVEC(9) EIGINV(3,2)=EIGVEC(9)*EIGVEC(4)-EIGVEC(7)*EIGVEC(6) EIGINV(3,3)=EIGVEC(7)*EIGVEC(5)-EIGVEC(8)*EIGVEC(4) WW=EIGVEC(7)*EIGINV(1,1)+EIGVEC(4)*EIGINV(2,1) 1 +EIGVEC(1)*EIGINV(3,1) DO 3060 N=1,3 DO 3060 NN=1,3 EIGINV(N,NN)=EIGINV(N,NN)/WW 3060 CONTINUE c DO 1000 N=1,NATMS CORDX(N)=CORD(N,1) CORDY(N)=CORD(N,2) CORDZ(N)=CORD(N,3) 1000 CONTINUE c DO 997 N=1,NATMS CARTDF(N,1)=CORDX(N)*EIGINV(1,1)+CORDY(N)*EIGINV(1,2)+ 1 CORDZ(N)*EIGINV(1,3) CARTDF(N,2)=CORDX(N)*EIGINV(2,1)+CORDY(N)*EIGINV(2,2)+ 1 CORDZ(N)*EIGINV(2,3) CARTDF(N,3)=CORDX(N)*EIGINV(3,1)+CORDY(N)*EIGINV(3,2)+ 1 CORDZ(N)*EIGINV(3,3) 997 CONTINUE 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_____________________________________________________________________________