$debug c------------------------------------------------------------------------------ c c GAMDMA - Conversion of DMA output from GAMESS into form suitable for C MIN16 and ALL c c NOTE: 1/ GAMESS prints out distributed moments NOT multipoles, hence c second and third moments have to be converted into c quadrupoles and octopoles c 2/ units in the STONE output section of GAMESS are all atomic units, c as against the conventional units in the ELMOM output c 3/ GAMDMA will convert the first set of STONE results it finds, co c take care when you want the MP2 and not SCF multipoles c 4/ GAMDMA locates the STONE output by finding the 'NET CHARGES' line c so only this part of the .OUT file needs to be available c c------------------------------------------------------------------------------ C C IMPORTANT: C C It is always a good idea to check the accuracy of conversion by recovering C the molecular multipoles with a program like MIN16 and comparing with those C in the GAMESS output. However: C C - total moments in the GAMESS output are relative to the origin at C coordinates labelled "POINT 1 X Y Z (BOHR) CHARGE" which is C where GAMESS thinks the centre of mass is C - DMA's from GAMDMA and thus resulting molecular multipoles are relative C to X,Y,X=0,0,0 C C Hence if numerical values in the "ELECTROSTATIC MOMENTS" output in GAMESS C are to be reproduced and the expansion point in GAMESS is not the C Cartesian origin, then a shift of the DMA's out of GAMDMA has to be made. C C This is done with GAMDMA -> "simple translation" option. Specify DeltaX, C DeltaY,DeltaZ as -X,-Y,-Z from the "POINT 1 X Y Z (BOHR) CHARGE" C line, and remember that X,Y,Z listed in GAMESS are in BOHR, while GAMDMA C wants DeltaX etc. in Angstroms. C C To obtain a DMA with cartesians as in GAMESS and recovering centre of mass C molecular multipoles (also as in GAMESS) it is necessary to do with DMAROT C first "simple shift of Cartesians" by -X,-Y,-Z followed C by "shift of multipole origin" by X,Y,Z c c C ver. 12.IV.2003 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c C Modification history: C C 10.12.97: creation C 31.03.00: modifications for PC-GAMESS C 28.11.01: modification for expansions other than with the ATOMS keyword C 12.04.03: elimination of bug in charge conversion and more commenting C______________________________________________________________________________ c c implicit real*8 (a-h,o-z) parameter (maxat=50) c character filnam*30,line*155 character*9 sitnam(maxat),atnam(maxat),ver real*8 itreat(maxat),del(3,3) common /blk1/x(maxat,3),xa(maxat,3),q(maxat),qa(maxat), * d(maxat,3),t(maxat,3,3),u(maxat,3,3,3) common /blk2/sitnam,filnam ver='OLDGAMESS' c C natoms - number of atoms in the molecule C nsites - number of DMA expansion sites C xa, qa, atnam - coordinates, charges and names of atoms C x, q, sitnam - coordinates, charges and names of DMA expansion sites C d, t, u - cartersian dipoles, quadrupoles, octopoles for DMA sites c write (*,297) 297 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | GAMDMA - Conversion of STONE group output from GAMESS', * ' into spherical DMA',T79,'|'/ * ' | components suitable for the electrostatic model', * T79,'|'/ * ' |',76(1H_),'|'/' version 12.IV.2003',T64,'Zbigniew KISIEL'/) c 1 write(*,'(1x//'' GAMESS file name: '',$)') read(*,'(a)',err=1)filnam open(2,file=filnam,status='old',err=1) c C...Locate STONE output C 2 read(2,'(A)',ERR=3,END=3)LINE IF(LINE(1:16).NE.' NET CHARGES')GOTO 2 C READ(2,'(A)')LINE READ(2,'(A)')LINE c c...Read information on atoms c natoms=0 44 READ(2,'(A)')LINE if(line(9:9).eq.'Z')then natoms=natoms+1 atnam(natoms)(1:9)=line(9:17) read(line(19:64),7) * qa(natoms),xa(natoms,1),xa(natoms,2),xa(natoms,3) 7 format(f9.5,f13.5,f12.5,f12.5) itreat(natoms)=1 goto 44 endif write(*,'(1x/'' Number of atoms:'',i5)')natoms c C...Identify DMA sites, their charges and coordinates c nsites=1 read(line(19:64),7) * q(nsites),x(nsites,1),x(nsites,2),x(nsites,3) read(line(9:17),'(a9)')sitnam(nsites) c 4 READ(2,'(A)')LINE if(line(22:22).eq.'.')then nsites=nsites+1 read(line(19:64),7) * q(nsites),x(nsites,1),x(nsites,2),x(nsites,3) read(line(9:17),'(a9)')sitnam(nsites) goto 4 endif write(*,'( '' Number of DMA sites:'',i5)')nsites c C...Check which sites are on atoms and adjust charge accordingly c do 45 i=1,nsites do 46 j=1,natoms if(x(i,1).eq.xa(j,1).and.x(i,2).eq.xa(j,2) . .and.x(i,3).eq.xa(j,3))then q(i)=q(i)+qa(j) itreat(j)=0 endif 46 continue 45 continue C c...dipoles c 5 READ(2,'(A)')LINE READ(2,'(A)')LINE READ(2,'(A)')LINE if(line(4:4).eq.'N')ver='PC-GAMESS' write(*,'(1x/'' GAMESS version is: '',a)')ver c n=0 8 read(2,'(a)')line n=n+1 if(ver.eq.'OLDGAMESS')then read(line(13:45),9 )d(n,1),d(n,2),d(n,3) else read(line(11:52),29)d(n,1),d(n,2),d(n,3) endif 9 format(6f11.5) 29 format(20f14.5) c 21 if(n.lt.nsites)goto 8 C c...quadrupoles: the GAMESS output in both .OUT and .PUN files is cartesian c second moments and a first step is to convert to cartesian quadrupoles c READ(2,'(A)')LINE READ(2,'(A)')LINE READ(2,'(A)')LINE READ(2,'(A)')LINE n=0 10 read(2,'(a)')line n=n+1 if(ver.eq.'OLDGAMESS')then read(line(13:78),9)t(n,1,1),t(n,2,2),t(n,3,3), * t(n,1,2),t(n,1,3),t(n,2,3) else read(line(11:94),29)t(n,1,1),t(n,2,2),t(n,3,3), * t(n,1,2),t(n,1,3),t(n,2,3) endif C qxx=0.5d0*(2.d0*t(n,1,1)-t(n,2,2)-t(n,3,3)) qyy=0.5d0*(2.d0*t(n,2,2)-t(n,1,1)-t(n,3,3)) qzz=0.5d0*(2.d0*t(n,3,3)-t(n,1,1)-t(n,2,2)) t(n,1,1)=qxx t(n,2,2)=qyy t(n,3,3)=qzz t(n,1,2)=1.5d0*t(n,1,2) t(n,1,3)=1.5d0*t(n,1,3) t(n,2,3)=1.5d0*t(n,2,3) c if(n.lt.nsites)goto 10 C c...octopoles: the GAMESS output in both .OUT and .PUN files is cartesian c third moments and a first step is to convert to cartesian octopoles c READ(2,'(A)')LINE READ(2,'(A)')LINE READ(2,'(A)')LINE READ(2,'(A)')LINE if(ver.eq.'OLDGAMESS')READ(2,'(A)')LINE n=0 12 read(2,'(a)')line n=n+1 if(ver.eq.'OLDGAMES')then read(line(13:78),9,err=3)u(n,1,1,1),u(n,2,2,2),u(n,3,3,3), * u(n,1,1,2),u(n,1,1,3),u(n,1,2,2) read(2,'(a)')line read(line(13:56),9,err=3)u(n,2,2,3),u(n,1,3,3),u(n,2,3,3), * u(n,1,2,3) else read(line(11:150),29,err=3)u(n,1,1,1),u(n,2,2,2),u(n,3,3,3), * u(n,1,1,2),u(n,1,1,3),u(n,1,2,2), * u(n,2,2,3),u(n,1,3,3),u(n,2,3,3), * u(n,1,2,3) endif c txxx=u(n,1,1,1) tyyy=u(n,2,2,2) tzzz=u(n,3,3,3) txxy=u(n,1,1,2) txxz=u(n,1,1,3) txyy=u(n,1,2,2) tyyz=u(n,2,2,3) txzz=u(n,1,3,3) tyzz=u(n,2,3,3) c ozzz=0.5d0*(5.d0*tzzz-3.d0*(tzzz+txxz+tyyz)) oxxz=0.5d0*(5.d0*txxz- (tzzz+txxz+tyyz)) oyyz=0.5d0*(5.d0*tyyz- (tzzz+txxz+tyyz)) oxxx=0.5d0*(5.d0*txxx-3.d0*(txxx+txyy+txzz)) oxyy=0.5d0*(5.d0*txyy- (txxx+txyy+txzz)) oxzz=0.5d0*(5.d0*txzz- (txxx+txyy+txzz)) oyyy=0.5d0*(5.d0*tyyy-3.d0*(tyyy+txxy+tyzz)) oxxy=0.5d0*(5.d0*txxy- (tyyy+txxy+tyzz)) oyzz=0.5d0*(5.d0*tyzz- (tyyy+txxy+tyzz)) c c...test code for printing cartesian octopoles c write(*,25)' n = ',n,' zzz',ozzz,' xxz',oxxz,' yyz',oyyz, c * ' xxx',oxxx,' xyy',oxyy,' xzz',oxzz, c * ' yyy',oyyy,' xxy',oxxy,' yzz',oyzz c25 format(1x,a5,i2,5(a,f8.4)/8x,5(a,f8.4)) c u(n,1,1,1)=oxxx u(n,2,2,2)=oyyy u(n,3,3,3)=ozzz u(n,1,1,2)=oxxy u(n,1,1,3)=oxxz u(n,1,2,2)=oxyy u(n,2,2,3)=oyyz u(n,1,3,3)=oxzz u(n,2,3,3)=oyzz u(n,1,2,3)=2.5d0*u(n,1,2,3) c if(n.lt.nsites)goto 12 c close(2) c c C...Deal with atoms, which are not used as expansion sites - subtract C contribution from the charge of such atom from DMA paramaters of C the nearest site c c do 98 i=1,3 do 99 j=1,3 del(i,j)=0.0d0 99 continue del(i,i)=1.0d0 98 continue c do 50 j=1,natoms if(itreat(j).eq.1)then rr=100.d0 do 51 i=1,nsites r=dsqrt( (xa(j,1)-x(i,1))**2 + (xa(j,2)-x(i,2))**2 . + (xa(j,3)-x(i,3))**2 ) if(r.lt.rr)then rr=r jr=i endif 51 continue c c...Corrections at site JR for charge on atom J c c...correct the charges q(jr)=q(jr)+qa(j) c c...correct the dipoles do 52 k=1,3 d(jr,k)=d(jr,k)+qa(j)*(xa(j,k)-x(jr,k)) 52 continue c c...correct the quadrupoles do 53 k=1,3 do 53 kk=1,3 t(jr,k,kk)=t(jr,k,kk) . +1.5d0*qa(j)*(xa(j,k)- x(jr,k)) . *(xa(j,kk)-x(jr,kk)) . -0.5d0*del(k,kk)*qa(j)*rr**2 53 continue c c...correct the octopoles do 54 k=1,3 do 54 kk=1,3 do 54 kkk=1,3 term=(xa(j,k) - x(jr,k)) *del(kk,kkk) . +(xa(j,kk) -x(jr,kk)) *del(k,kkk) . +(xa(j,kkk)-x(jr,kkk)) *del(k,kk) term3=-0.5d0*qa(j)*(rr**2*term . -5.0d0*(xa(j,k) -x(jr,k)) . *(xa(j,kk) -x(jr,kk)) . *(xa(j,kkk)-x(jr,kkk)) ) u(jr,k,kk,kkk)=u(jr,k,kk,kkk)+term3 54 continue c itreat(j)=0 endif 50 continue c c...convert to sphericals and output c call consph(nsites) stop c 3 WRITE(*,'(1X//'' ****ERROR reading input:''/1x,a/)')line C stop end c c---------------------------------------------------------------------------- c SUBROUTINE CONSPH(nsites) c c Conversion of cartesian multipoles into spherical ones and c on the fly file output. Sphericals are in standard order: c c q00 c q10 q11c q11s c q20 q21c q21s q22c q22s c q30 q31c q31s q32c q32s q33c q33s c q40 q41c q41s q42c q42s q43c q43s q44c q44s c c Indices for cartesians are assumed to be X,Y,Z for 1,2,3 resp. c implicit real*8 (a-h,o-z) parameter (maxat=50) c real*8 out(9) character sitnam(maxat)*9,filnam*30 common /blk1/x(maxat,3),xa(maxat,3),q(maxat),qa(maxat), * d(maxat,3),t(maxat,3,3),u(maxat,3,3,3) common /blk2/sitnam,filnam c c OPEN(3,FILE='gamdma.out',STATUS='unknown') c write(3,'(i4,'' DMA from file '',a)')nsites,filnam c do 1 n=1,nsites c c...coordinates c write(3,207)x(n,1),x(n,2),x(n,3),sitnam(n),n 207 format(3f10.5,10x,a9,'<- DMA site',i3) c c...charge c write(3,206)q(n) c c...dipole c write(3,206)d(n,3),d(n,1),d(n,2) c c...quadrupole c out(1)=t(n,3,3) out(2)=dsqrt(4.d0/3.d0)*t(n,1,3) out(3)=dsqrt(4.d0/3.d0)*t(n,2,3) out(4)=dsqrt(1.d0/3.d0)*(t(n,1,1)-t(n,2,2)) out(5)=dsqrt(4.d0/3.d0)*t(n,1,2) write(3,206)(out(i),i=1,5) c c...octopole c out(1)=u(n,3,3,3) out(2)=dsqrt(3.d0/2.d0)*u(n,1,3,3) out(3)=dsqrt(3.d0/2.d0)*u(n,2,3,3) out(4)=dsqrt(3.d0/5.d0)*(u(n,1,1,3)-u(n,2,2,3)) out(5)=2.d0*dsqrt(3.d0/5.d0)*u(n,1,2,3) out(6)=dsqrt(1.d0/10.d0)*(u(n,1,1,1)-3.d0*u(n,1,2,2)) out(7)=dsqrt(1.d0/10.d0)*(3.d0*u(n,1,1,2)-u(n,2,2,2)) write(3,206)(out(i),i=1,7) c c...hexadecapole c write(3,'(5x,a)')'0,0,0,0,0,0,0,0,0' c 206 format(9f10.5) c 1 continue c close(3) write(*,'(1x//'' -----Output written to GAMDMA.OUT''//)') c return end c c---------------------------------------------------------------------------- c----------------------------------------------------------------------------