Electronic Supplementary Material for CrystEngComm This journal is © The Royal Society of Chemistry 2003 c c PIXEL.II Program Package may 2003 version c c This is the main module for the calculation of energies c c THIS IS FREE SOFTWARE FOR ACADEMIC AND NON-PROFIT ORGANIZATIONS c USE AND REDISTRIBUTION OF THIS SOFTWARE c FOR PROFIT PURPOSES IS FORBIDDEN c program pixel2 c c A.Gavezzotti c Dipartimento di Chimica Strutturale e Stereochimica Inorganica c University of Milano c via Venezian 21 c I-20133 Milano (Italy) c angelo.gavezzotti at unimi.it c c uses a GAUSSIAN-98 electron density file (binary format) c for a box containing nmsolu solute molecules c and nmsolv solvent molecules c calculates electrostatic energy, electrostatic field, c polarization energy, overlap repulsion, c and dispersion energy c c charges in electrons, coordinates in angstrom, c polarizabilities in cubic angstrom c ionization potential in a.u. c implicit real*8 (a-h,o-z) character*4 titl(10) character*2 flab c last letter = u : solute quantities c last letter = V : solvent quantities dimension epse(60000,3),bufe(60000,3) dimension epstot(3),debuf(500),boxbuf(500),polze(30) c c xze is the starting point for the density grid in angstrom c fze is the atomic number (less core electrons if valence density) c elto is the nuclear charge c x are atomic nuclear coordinates in angstrom,q atomic charges c atomic species indicators (see manual) dimension xzeu(3),xzev(3),fzeu(50),fzev(50), 1 xu(50,3),xv(50,3),qu(50),qv(50),ispu(50),ispv(50) c for each atomic species, ravdw are van der waals atomic radii, c zval and ztot are atomic numbers (valence and total) dimension ravdw(30),zval(30),ztot(30) c polat is atomic polarizability dimension polatu(50),polatv(50) c ro is the main electronic property table c this table has x,y,z (angstrom), ro (electrons), c polarizability,(10-30 cubic angstrom),original density dimension rou(60000,6),rov(60000,6) c npat is number of ro points for each atomic basin dimension npatu(50),npatv(50) c grid tables have the center of mass position and rotation angles c for each molecule in the box dimension gridu(500,6),gridv(500,6) c tipo is the 6-exp parameter table COMMON /A/ tipo(30,30,3),rinter(30),rintm(30,30),flab(30) common /b/ eltot,eunir,eunia,eqq,edisp,erep,icry c tcry and fcry are matrices and vectors for cluster c molecules in the case of a crystal common /cry/ fcry(500,3,3),tcry(500,3) common /ovla/ bdrep,bdexp,collis,ddamp c box arrays contain the complete ro density map common /boxes/ boxu(5000000),boxv(5000000), 1 boxmiu(3),boxmau(3),boxmiv(3),boxmav(3), 2 rominu(3),romaxu(3),rominv(3),romaxv(3) c c library atomic radii c H C N O F Si P S Cl ar - DATA RAVDW/1.10,0.0,1.77,1.64,1.58,1.46,2.0,1.9,1.81,1.76,1.8,0.0, c Br - I O O O N N S O H H H H 1 1.87,0.0,2.03,3*1.58,1.64,1.64,1.81,1.58,2*0.0,4*1.10,0.,0./ c intermolecular contact radii c H ar C N O F Si P S Cl - - DATA rinter/1.1,2.7,1.77,1.64,1.58,1.46,2.0,1.9,1.81,1.76,1.8,0.0, c Br - I O O O N N S O H H H H 1 1.87,0.0,2.03,3*1.30,1.30,1.30,1.81,1.30,2*0.0,4*0.30,0.,0./ c valence and total atomic numbers data ztot/1.,0.,6.,7.,8.,9.,14.,15.,16.,17.,18.0,0.0, 1 35.,0.0,53.,3*8.0,7.0,7.0,16.,8.0,2*0.0,4*1.0,0.,0./ data zval/1.,8.,4.,5.,6.,7.,4.,5.,6.,7.,8.,0., 1 7.,0.,7.,3*6.,5.,5.,6.,6.,0.,0.,4*1.,0.,0./ c c library atomic polarizabilities c H C N O F Si - S Cl ar - data polze/0.39,0.0,1.05,0.88,0.80,0.41,0.0,0.0,3.00,2.30,1.7,0.0, c Br - I O O O N N S O H(O) H(N) 1 3.27,0.0,5.50,3*0.80,2*0.88,3.7,0.80,2*0.0,2*0.39,2*0.39,2*0.0/ data flab/'H ','XX','C ','N ','O ','F ','SI','P ','S ','CL', 1 'AR','ZZ','BR','WW','I ','O ','O ','O ','N ','N ','S ', 2 'O ',' ',' ','H ','H ','H ','H ',' ',' '/ c open files c unit 16 output file with molecule-molecule energies c unit 8 printout c unit 9 card input open (unit=16, file='pixel.mlc') open (unit=9, file='pixel.inp') open (unit=8, file='pixel.oxp') c unit 12 scratch for fields and potentials c unit 14 scratch file with condensed ro open (unit=14, file='fort.14') open (unit=12, form='unformatted',file='fort.12') c unit 21 scratch in generation of condensed charge file open (unit=21, form='unformatted',file='fort.21') c state constants c nro max charge points, nre max grid steps in density file c maxat max no. of atoms in molecules, solvent or solute c maxmol max no. of molecules in box, solvent or solute c nreb max grid steps in one dimension nro=60000 nre=5000000 nreb=500 maxat=50 maxmol=500 c compute conversion factors avoga=6.0221367d+23 bohr=0.52917d+00 coule=1.6021773d-19 pigr=3.1415926536d+00 pz=2*pigr/360.0 epszer=8.85419d-12 angmet=1.0d+10 fourp=4.0*pigr*epszer facpot=coule*angmet/fourp facen=facpot*coule facfie=facpot*angmet read(9,791) titl 791 format(1x,10a4) write(6,790) titl write(8,790) titl 790 format(//2x,'PIXEL.II, may 2003 version'/2x,10a4) c c general card input for the whole run c read 'collision parameter', damping factor for e(disp) c and damping parameter for e(pol) c bdrep and bdexp are two scaling parameters for repulsion read(9,*) collis,ddamp,fiemax,bdrep,bdexp c sets defaults for parameters c default for collis will be set after steps are read rid=0.000001 c these are the four parameters of the theory if(bdrep.lt.rid) bdrep=3700. if(ddamp.lt.rid) ddamp=3.25 if(bdexp.lt.rid) bdexp=1.00 if(fiemax.lt.rid) fiemax=200. write(8,788) ddamp,fiemax,bdrep,bdexp write(6,788) ddamp,fiemax,bdrep,bdexp 788 format(/' == Scaling, Edisp,Epol,Erep(tot-exp)'/12x,f8.3, 1 f9.2,f9.2,f7.3) c c read screenout parameters c read condensation level for solute c ival=0 valence density, ne.0 total density c then min and max renormalized charge read(9,*) nrdu,ivalu,romiu,romau if(nrdu.eq.0) nrdu=4 if(romiu.eq.0.0) romiu=1.0d-06 if(romau.lt.rid) romau=9999. if(nrdu.le.10) go to 10 write(6,212) nrdu 212 format(2x,'ro reduction number too high',i5) stop 10 continue c screenout parameters for solvent read(9,*) nrdv,ivalv,romiv,romav if(nrdv.eq.0) nrdv=4 if(romiv.eq.0.0) romiv=1.0d-06 if(romav.lt.rid) romav=9999. if(nrdv.le.10) go to 1111 write(6,212) nrdv stop 1111 continue rewind 12 rewind 14 rewind 21 c c now read card input proper of one molecular system c number of solute and solvent molecules c number of atoms in solute and solvent molecule c nmsolu must be >0, nmsolv=0 is ok c if nmsolv<0 generate a crystal structure read(9,*) nmsolu,nmsolv, nasolu,nasolv c load 6-exp potential energy parameters call uni6 c set icry=0 if nocrystal, = 1 if crystal icry=0 if(nmsolv.lt.0) icry=1 if(nmsolu.le.maxmol.and.nmsolv.le.maxmol) go to 7 write(6,8) nmsolu,nmsolv 8 format(2x,'emergency stop',2i5) stop 7 if(nasolu.le.maxat.and.nasolv.le.maxat) go to 9 write(6,8) nasolu,nasolv stop 9 continue c read average ionization potential and charge for solute read(9,*) avionu,chargu if(avionu.lt.rid) avionu=1.0 write(8,849) avionu,chargu 849 format(' == Solute ion. potential and charge',f10.4,f8.2) c c start reading electron density file for solute open (unit=11, form='unformatted', file='solu.bin') c c reads a Gaussian electron density for solute molecule c c starting point for electron density (bohr) read(11) idum,(xzeu(i),i=1,3) c603 format(i5,4f12.0) c number of steps, step (bohr) x, y, z read(11) nxu,dxxu,dxyu,dxzu read(11) nyu,dyxu,dyyu,dyzu read(11) nzu,dzxu,dzyu,dzzu nxuor=nxu nyuor=nyu nzuor=nzu dxxu=dxxu*bohr dyyu=dyyu*bohr dzzu=dzzu*bohr xzeu(1)=xzeu(1)*bohr xzeu(2)=xzeu(2)*bohr xzeu(3)=xzeu(3)*bohr stpu=(dxxu+dyyu+dzzu)/3.0 ntox=nxu*nyu*nzu if(nxu.gt.nreb.or.nyu.gt.nreb.or.nzu.gt.nreb) go to 6334 if(ntox.le.nre) go to 634 6334 write(6,633) nxu,nyu,nzu 633 format(2x,'emergency stop, too many density pts',3i6) stop 634 continue boxmiu(1)=xzeu(1) boxmiu(2)=xzeu(2) boxmiu(3)=xzeu(3) boxmau(1)=boxmiu(1)+dfloat(nxu-1)*dxxu boxmau(2)=boxmiu(2)+dfloat(nyu-1)*dyyu boxmau(3)=boxmiu(3)+dfloat(nzu-1)*dzzu c read atomic coordinates (bohr) of solute in ed reference frame write(8,639) 639 format(/' == Solute atoms Z x,y,z charge polar. vdW R') fnucu=0.0 ridic=0.0001 do 604 i=1,nasolu read(11) idum,edum,(xu(i,k),k=1,3) c read atomic species, charge, atomic polarizability read(9,*) nax,ispen,qu(i),polatu(i) 663 format(i5,4f12.0) if(ivalu.ne.0) fzeu(i)=ztot(ispen) if(ivalu.eq.0) fzeu(i)=zval(ispen) ispu(i)=ispen fnucu=fnucu+fzeu(i) if(polatu(i).lt.ridic) polatu(i)=polze(ispen) do 605 k=1,3 605 xu(i,k)=xu(i,k)*bohr write(8,631) ispen,fzeu(i),(xu(i,k),k=1,3),qu(i),polatu(i), 1 ravdw(ispen) 631 format(2x,i4,1x,f7.2,2x,3f10.5,f8.4,f7.2,1x,f7.2) 604 continue c first read uncondensed ro matrix bma=-99999. bmi=99999. do 6671 i=1,nxu do 6671 j=1,nyu read(11) (boxbuf(k),k=1,nzu) do 6672 k=1,nzu ibox=(i-1)*nyu*nzu+(j-1)*nzu+k boxu(ibox)=boxbuf(k)/bohr**3 if(boxu(ibox).gt.bma) bma=boxu(ibox) if(boxu(ibox).lt.bmi) bmi=boxu(ibox) 6672 continue 6671 continue c write(8,6673) bmi,bma 6673 format(2x,'solute,min and max original density',2e12.4) c then condense call sumden (boxu,nxu,nyu,nzu,nrdu,dxxu,dyyu,dzzu) c electron density points are now charge points totel=0.0 irou=0 totelu=0.0 debmen=0.0 debplu=0.0 xzero=xzeu(1)+0.5*dfloat(nrdu-1)*dxxu stepx=dfloat(nrdu)*dxxu yzero=xzeu(2)+0.5*dfloat(nrdu-1)*dyyu stepy=dfloat(nrdu)*dyyu zzero=xzeu(3)+0.5*dfloat(nrdu-1)*dzzu stepz=dfloat(nrdu)*dzzu stepu=(stepx+stepy+stepz)/3.0 istop=0 do 695 i=1,3 rominu(i)=999999. 695 romaxu(i)=-999999. do 606 i=1,nxu do 607 j=1,nyu read(14,608) (debuf(k),k=1,nzu) 608 format(6e13.5) vcelu=stepx*stepy*stepz c screenout do 609 k=1,nzu deb=debuf(k) totel=totel+deb if((deb.gt.romiu).and.(deb.lt.romau)) go to 11 if(deb.lt.romiu) debmen=debmen+deb if(deb.gt.romau) debplu=debplu+deb go to 609 11 totelu=totelu+deb irou=irou+1 if(irou.le.nro) go to 12 istop=1 go to 609 12 rou(irou,1)=xzero+dfloat(i-1)*stepx if(rou(irou,1).lt.rominu(1)) rominu(1)=rou(irou,1) if(rou(irou,1).gt.romaxu(1)) romaxu(1)=rou(irou,1) rou(irou,2)=yzero+dfloat(j-1)*stepy if(rou(irou,2).lt.rominu(2)) rominu(2)=rou(irou,2) if(rou(irou,2).gt.romaxu(2)) romaxu(2)=rou(irou,2) rou(irou,3)=zzero+dfloat(k-1)*stepz if(rou(irou,3).lt.rominu(3)) rominu(3)=rou(irou,3) if(rou(irou,3).gt.romaxu(3)) romaxu(3)=rou(irou,3) c charge rou(irou,4)=-deb 609 continue 607 continue 606 continue write(6,611) nxu,nyu,nzu,stepx,stepy,stepz,vcelu write(8,611) nxu,nyu,nzu,stepx,stepy,stepz,vcelu 611 format(' == solutes,steps,vol, (A)',/2x,3i5,3f9.4,f10.5) write(8,612) totel,totelu,irou write(6,612) totel,totelu,irou 612 format(' == no. electrons',f10.5,2x,'remaining-pixels', 1 2x,f10.5,i7) write(8,693) romiu,romau,debmen,debplu write(6,693) romiu,romau,debmen,debplu 693 format(' == q min and max',e12.4,f10.2, 1 /' == electrons out low and high',2e12.4) c output full box limits c write(8,694) boxmiu(1),boxmau(1),stepu,boxmiu(2),boxmau(2), c 1 stepu,boxmiu(3),boxmau(3),stepu c694 format(1x,'full box limits in x,y,z, steps'/2x,3(2x,3f8.3)) c write(8,697) rominu(1),romaxu(1),rominu(2),romaxu(2), c 1 rominu(3),romaxu(3) 697 format(1x,'red. box limits in x,y,z'/2x,6f9.3) if(istop.ne.0) write(6,13) irou 13 format(2x,'emergency stop ',i8) if(istop.ne.0) stop c search atomic basin to which charge point belongs do 674 ii=1,nasolu 674 npatu(ii)=0 fnucuq=fnucu-chargu check=dabs(dabs(totel)-dabs(fnucuq)) ch=0.1*dabs(totel) if(check.lt.ch) go to 1776 write(6,1775) totel,fnucuq 1775 format(' something wrong with density',2e12.4) stop 1776 continue factor=fnucuq/totelu rototu=0.0 do 672 i=1,irou rpomi=9999999. do 711 kp=1,nasolu rpo=(rou(i,1)-xu(kp,1))**2 rpo=rpo+(rou(i,2)-xu(kp,2))**2 rpo=rpo+(rou(i,3)-xu(kp,3))**2 rpo=sqrt(rpo) ispek=ispu(kp) rats=rpo/ravdw(ispek) if(rats.gt.rpomi) go to 711 kmin=kp rpomi=rats 711 continue c factor rescales to exact number of electrons rou(i,4)=rou(i,4)*factor rototu=rototu+rou(i,4) npatu(kmin)=npatu(kmin)+1 polact=polatu(kmin) c polarizability rou(i,5)=-polact*rou(i,4)/fzeu(kmin) c original density rou(i,6)=-rou(i,4)/vcelu 672 continue write(6,676) fnucu,rototu write(8,676) fnucu,rototu 676 format(' == check renormalized total charges',2x,2f12.6) c write(8,675) (k,npatu(k),k=1,nasolu) c675 format(' == solute,no. q points per atom'/2x,(12(i3,i5))) avionv=avionu if(nmsolv.le.0) go to 701 c c reads a Gaussian electron density for solvent molecule c read average ionization potential and charge for solvent read(9,*) avionv,chargv if(avionv.lt.rid) avionv=1.0 write(8,850) avionv,chargv 850 format(' == Solvent ion. potential and charge',f10.4,f8.2) open (unit=15, form='unformatted', file='solv.bin') read(15) idum,(xzev(i),i=1,3) read(15) nxv,dxxv,dxyv,dxzv read(15) nyv,dyxv,dyyv,dyzv read(15) nzv,dzxv,dzyv,dzzv nxvor=nxv nyvor=nyv nzvor=nzv dxxv=dxxv*bohr dyyv=dyyv*bohr dzzv=dzzv*bohr stpv=(dxxv+dyyv+dzzv)/3.0 xzev(1)=xzev(1)*bohr xzev(2)=xzev(2)*bohr xzev(3)=xzev(3)*bohr ntox=nxv*nyv*nzv if(nxv.gt.nreb.or.nyv.gt.nreb.or.nzv.gt.nreb) go to 6335 if(ntox.le.nre) go to 635 6335 write(6,633) nxv,nyv,nzv stop 635 continue boxmiv(1)=xzev(1) boxmiv(2)=xzev(2) boxmiv(3)=xzev(3) boxmav(1)=boxmiv(1)+dfloat(nxv-1)*dxxv boxmav(2)=boxmiv(2)+dfloat(nyv-1)*dyyv boxmav(3)=boxmiv(3)+dfloat(nzv-1)*dzzv write(8,638) 638 format( ' == Solvent atoms') fnucv=0.0 ridic=0.0001 do 614 i=1,nasolv read(15) izev,elto,(xv(i,k),k=1,3) read(9,*) nax,ispen,qv(i),polatv(i) if(ivalv.ne.0) fzev(i)=ztot(ispen) if(ivalv.eq.0) fzev(i)=zval(ispen) ispv(i)=ispen fnucv=fnucv+fzev(i) if(polatv(i).lt.ridic) polatv(i)=polze(ispen) do 615 k=1,3 615 xv(i,k)=xv(i,k)*bohr write(8,631) ispen,fzev(i),(xv(i,k),k=1,3),qv(i),polatv(i), 1 ravdw(ispen) 614 continue bma=-99999. bmi=99999. do 6675 i=1,nxv do 6675 j=1,nyv read(15) (boxbuf(k),k=1,nzv) do 6677 k=1,nzv ibox=(i-1)*nyv*nzv+(j-1)*nzv+k boxv(ibox)=boxbuf(k)/bohr**3 if(boxv(ibox).gt.bma) bma=boxv(ibox) if(boxv(ibox).lt.bmi) bmi=boxv(ibox) 6677 continue 6675 continue c write(8,6678) bmi,bma c6678 format(2x,'solvent,min and max original density',2e12.4) call sumden (boxv,nxv,nyv,nzv,nrdv,dxxv,dyyv,dzzv) totel=0.0 irov=0 totelv=0.0 debmen=0.0 debplu=0.0 xzero=xzev(1)+0.5*dfloat(nrdv-1)*dxxv stepx=dfloat(nrdv)*dxxv yzero=xzev(2)+0.5*dfloat(nrdv-1)*dyyv stepy=dfloat(nrdv)*dyyv zzero=xzev(3)+0.5*dfloat(nrdv-1)*dzzv stepz=dfloat(nrdv)*dzzv stepv=(stepx+stepy+stepz)/3.0 istop=0 do 696 i=1,3 rominv(i)=999999. 696 romaxv(i)=-999999. do 616 i=1,nxv do 617 j=1,nyv read(14,608) (debuf(k),k=1,nzv) vcelv=stepx*stepy*stepz do 619 k=1,nzv deb=debuf(k) totel=totel+deb if((deb.gt.romiv).and.(deb.lt.romav)) go to 101 if(deb.lt.romiv) debmen=debmen+deb if(deb.gt.romav) debplu=debplu+deb go to 619 101 totelv=totelv+deb irov=irov+1 if(irov.le.nro) go to 102 istop=1 go to 619 102 rov(irov,1)=xzero+dfloat(i-1)*stepx if(rov(irov,1).lt.rominv(1)) rominv(1)=rov(irov,1) if(rov(irov,1).gt.romaxv(1)) romaxv(1)=rov(irov,1) rov(irov,2)=yzero+dfloat(j-1)*stepy if(rov(irov,2).lt.rominv(2)) rominv(2)=rov(irov,2) if(rov(irov,2).gt.romaxv(2)) romaxv(2)=rov(irov,2) rov(irov,3)=zzero+dfloat(k-1)*stepz if(rov(irov,3).lt.rominv(3)) rominv(3)=rov(irov,3) if(rov(irov,3).gt.romaxv(3)) romaxv(3)=rov(irov,3) rov(irov,4)=-deb 619 continue 617 continue 616 continue write(6,621) nxv,nyv,nzv,stepx,stepy,stepz,vcelv write(8,621) nxv,nyv,nzv,stepx,stepy,stepz,vcelv 621 format(' == solvents,steps,vol, (A)',/2x,3i5,3f9.4,f10.5) write(8,612) totel,totelv,irov write(6,612) totel,totelv,irov write(8,693) romiv,romav,debmen,debplu write(6,693) romiv,romav,debmen,debplu c output full box limits c write(8,694) boxmiv(1),boxmav(1),boxmiv(2),boxmav(2), c 1 boxmiv(3),boxmav(3) c write(8,697) rominv(1),romaxv(1),rominv(2),romaxv(2), c 1 rominv(3),romaxv(3) if(istop.ne.0) write(6,13) irov if(istop.ne.0) stop do 678 ii=1,nasolv 678 npatv(ii)=0 fnucvq=fnucv-chargv check=dabs(dabs(totel)-dabs(fnucvq)) ch=0.1*dabs(totel) if(check.lt.ch) go to 1777 write(6,1775) totel,fnucvq stop 1777 continue factor=fnucvq/totelv rototv=0.0 do 671 i=1,irov rpomi=9999999. do 712 kp=1,nasolv rpo=(rov(i,1)-xv(kp,1))**2 rpo=rpo+(rov(i,2)-xv(kp,2))**2 rpo=rpo+(rov(i,3)-xv(kp,3))**2 rpo=sqrt(rpo) ispek=ispv(kp) rats=rpo/ravdw(ispek) if(rats.gt.rpomi) go to 712 kmin=kp rpomi=rats 712 continue rov(i,4)=rov(i,4)*factor rototv=rototv+rov(i,4) npatv(kmin)=npatv(kmin)+1 polact=polatv(kmin) rov(i,5)=-polact*rov(i,4)/fzev(kmin) rov(i,6)=-rov(i,4)/vcelv 671 continue write(6,676) fnucv,rototv write(8,676) fnucv,rototv c write(8,680) (k,npatv(k),k=1,nasolv) 680 format(' == Solvent,no. q points per atom'/2x,(12(i3,i5))) 701 continue c c new system: new grid points for dimer case c or new crystal for the same molecule ridic=0.000001 if(stepv.lt.ridic) stepv=stepu stepav=0.5*(stepu+stepv) if(collis.lt.ridic) collis=stepav*0.5 write(6,789) collis,ddamp,fiemax write(8,789) collis,ddamp,fiemax 789 format(' ===== Start energy calculations ====='/, 1 2x,'collision parameter, damp factors for Edisp and Epol', 1 /2x,f10.4,1x,f10.3,1x,f10.1) 997 continue if(nmsolv.lt.0) nmsolu=0 write(16,1771) nmsolu,nmsolv 1771 format(2i5) rewind 12 if (icry.ne.0) go to 305 c non-crystal case write(8,641) 641 format(/,' Solute positions,c.o.m. and euler angles') c open interface file for schakal (E.Keller) open (unit=27, file='pixel.dat') do 301 i=1,nmsolu read(9,*) (gridu(i,k),k=1,6) write(8,331) i,(gridu(i,k),k=1,6) write(6,331) i,(gridu(i,k),k=1,6) 331 format(1x,'Grid vector mol.',i4,3f8.3,3f8.1,4x) 301 continue if(nmsolv.le.0) go to 303 c write(8,642) c642 format(' Solvent positions,c.o.m. and euler angles') do 304 i=1,nmsolv read(9,*) (gridv(i,k),k=1,6) write(8,331) i,(gridv(i,k),k=1,6) write(6,331) i,(gridv(i,k),k=1,6) 304 continue go to 303 305 continue c crystal, prepare space group matrices c read cutoff shell c first two numbers define a shell (min-max) c last number is top cutoff expected in a calculation read(9,*) cutmi,cutma,cuttop write(6,466) cutmi,cutma,cuttop write(8,466) cutmi,cutma,cuttop 466 format(/2x,'CRYSTAL; between',2f7.2,2x,'top cutoff',f7.2) call precry (cutmi,cutma,cuttop,nmsolu) write(6,467) nmsolu write(8,467) nmsolu 467 format(2x,i7,' molecules in crystal shell') 303 continue c c start energy calculations c call doall(1,nmsolu,nmsolu,nasolu,nasolu,xu,gridu,xu,gridu, 1 fzeu,fzeu,irou,rou,irou,rou,0,0,ispu,ispu,qu,qu, 2 avionu,avionu,boxu,boxu,boxmiu,boxmiu,boxmau,boxmau,stpu,stpu, 3 nyuor,nzuor,nyuor,nzuor) c compute total energies c ff force field energies c ec ed er ecr are coul, disp, rep, pixel energies ffreuu=eunir ffatuu=eunia fftouu=eunir+eunia ffqquu=eqq*angmet*coule*coule/fourp ffqquu=ffqquu*avoga*0.001 ecuu=eltot*avoga*0.001 eduu=edisp eruu=erep write(6,330) ffreuu,ffatuu,fftouu,ffqquu,ecuu,eduu,eruu c write(8,330) ffreuu,ffatuu,fftouu,ffqquu,ecuu,eduu,eruu 330 format(2x,'e(solute-solute):6EX,qiqj,coul,disp,rep', 1 /2x,8f8.1) c solute-solvent part c zero all solvent-related energies ffreuv=0.0 ffatuv=0.0 fftouv=0.0 ffqquv=0.0 ecuv=0.0 eduv=0.0 eruv=0.0 ffrevu=0.0 ffatvu=0.0 fftovu=0.0 ffqqvu=0.0 ecvu=0.0 edvu=0.0 ervu=0.0 ffrevv=0.0 ffatvv=0.0 fftovv=0.0 ffqqvv=0.0 ecvv=0.0 edvv=0.0 ervv=0.0 if(nmsolv.le.0) go to 335 call doall(2,nmsolu,nmsolv,nasolu,nasolv,xu,gridu,xv,gridv, 1 fzeu,fzev,irou,rou,irov,rov,0,1,ispu,ispv,qu,qv, 2 avionu,avionv,boxu,boxv,boxmiu,boxmiv,boxmau,boxmav,stpu,stpv, 3 nyuor,nzuor,nyvor,nzvor) ffreuv=eunir ffatuv=eunia fftouv=eunir+eunia ffqquv=eqq*angmet*coule*coule/fourp ffqquv=ffqquv*avoga*0.001 ecuv=eltot*avoga*0.001 eduv=edisp eruv=erep write(6,332) ffreuv,ffatuv,fftouv,ffqquv,ecuv,eduv,eruv c write(8,332) ffreuv,ffatuv,fftouv,ffqquv,ecuv,eduv,eruv 332 format(2x,'e(solute-solvnt):6EX,qiqj,coul,disp,rep', 1 /2x,4f8.1,f11.1,3f8.1) c solvent-solute part call doall(3,nmsolv,nmsolu,nasolv,nasolu,xv,gridv,xu,gridu, 1 fzev,fzeu,irov,rov,irou,rou,1,0,ispv,ispu,qv,qu, 2 avionv,avionu,boxv,boxu,boxmiv,boxmiu,boxmav,boxmau,stpv,stpu, 3 nyvor,nzvor,nyuor,nzuor) ffrevu=eunir ffatvu=eunia fftovu=eunir+eunia ffqqvu=eqq*angmet*coule*coule/fourp ffqqvu=ffqqvu*avoga*0.001 ecvu=eltot*avoga*0.001 edvu=edisp ervu=erep write(6,333) ffrevu,ffatvu,fftovu,ffqqvu,ecvu,edvu,ervu c write(8,333) ffrevu,ffatvu,fftovu,ffqqvu,ecvu,edvu,ervu 333 format(2x,'e(solvnt-solute):6EX,qiqj,coul,disp,rep', 1 /2x,4f8.1,f11.1,3f8.1) c solvent-solvent part call doall(1,nmsolv,nmsolv,nasolv,nasolv,xv,gridv,xv,gridv, 1 fzev,fzev,irov,rov,irov,rov,1,1,ispv,ispv,qv,qv, 2 avionv,avionv,boxv,boxv,boxmiv,boxmiv,boxmav,boxmav,stpv,stpv, 3 nyvor,nzvor,nyvor,nzvor) ffrevv=eunir ffatvv=eunia fftovv=eunir+eunia ffqqvv=eqq*angmet*coule*coule/fourp ffqqvv=ffqqvv*avoga*0.001 ecvv=eltot*avoga*0.001 edvv=edisp ervv=erep write(6,334) ffrevv,ffatvv,fftovv,ffqqvv,ecvv,edvv,ervv c write(8,334) ffrevv,ffatvv,fftovv,ffqqvv,ecvv,edvv,ervv 334 format(2x,'e(solvnt-solvnt):6EX,qiqj,coul,disp,rep', 1 /2x,8f8.1) 335 continue ffreto=ffreuu+ffreuv+ffrevu+ffrevv ffatto=ffatuu+ffatuv+ffatvu+ffatvv fftoto=fftouu+fftouv+fftovu+fftovv ffqqto=ffqquu+ffqquv+ffqqvu+ffqqvv ectot=ecuu+ecuv+ecvu+ecvv edtot=eduu+eduv+edvu+edvv ertot=eruu+eruv+ervu+ervv c c compute polarization energy c first sum up all electric field contributions c data on tape 12 are for molecules i-j c i-th molecule polarizing j-th molecule c c solute molecules uitota=0.0 nup=nmsolu if(icry.ne.0) nup=1 do 121 k=1,nup rewind 12 do 122 i=1,nro do 122 j=1,3 epstot(j)=0.0 122 epse(i,j)=0.0 103 read(12,end=110) i,iss1,j,iss2,nr2,elee,elen,elnn read(12) ((bufe(m,n),n=1,3),m=1,nr2) 321 format(4i4,i6,4e13.5) 323 format(5d16.8) c check that it is the k-th solute molecule if(k.ne.j) go to 103 if(iss2.ne.0) go to 103 do 104 m=1,irou do 104 n=1,3 epstot(n)=epstot(n)+bufe(m,n) 104 epse(m,n)=epse(m,n)+bufe(m,n) go to 103 c finished reading tape, electric field completed 110 continue uinda=0.0 unod=0.0 do 111 m=1,irou epsmod=epse(m,1)**2+epse(m,2)**2+epse(m,3)**2 epsi=dsqrt(epsmod) c compute polarization energy contribution in kJ/mol uparz=rou(m,5)*epsmod uparz=-2.0*pigr*epszer*0.001*avoga*1.0d-30*uparz unod=unod+uparz epshow=epsi*1.0d-10 epsh=epshow if(epsh.ge.fiemax) epsh=0.9999*fiemax damp=epsh/(fiemax-epsh) damp=dexp(-damp) epsdam=epsi*damp uparz=rou(m,5)*epsdam**2 uparz=-2.0*pigr*epszer*0.001*avoga*1.0d-30*uparz uinda=uinda+uparz c write(26,868) epshow,uparz,rou(m,5) 868 format(2f10.3,e13.4) 111 continue uitota=uitota+uinda write(6,112) k,uinda,unod write(8,112) k,uinda,unod 112 format(2x,'polariz.mol.',i4,' Ep,no-damp,at-field',3f8.2) 121 continue c c solvent molecules if(nmsolv.le.0) go to 999 c write(6,441) c441 format(2x,'computing solvent polarization') do 149 k=1,nmsolv rewind 12 do 148 i=1,nro do 148 j=1,3 epstot(j)=0.0 148 epse(i,j)=0.0 143 read(12,end=140) i,iss1,j,iss2,nr2,elee,elen,elnn read(12) ((bufe(m,n),n=1,3),m=1,nr2) if(k.ne.j) go to 143 if(iss2.ne.1) go to 143 do 144 m=1,irov do 144 n=1,3 epstot(n)=epstot(n)+bufe(m,n) 144 epse(m,n)=epse(m,n)+bufe(m,n) go to 143 140 continue uinda=0.0 unod=0.0 do 141 m=1,irov epsmod=epse(m,1)**2+epse(m,2)**2+epse(m,3)**2 uparz=rov(m,5)*epsmod uparz=-2.0*pigr*epszer*0.001*avoga*1.0d-30*uparz epshow=dsqrt(epsmod)*1.0d-10 epsh=epshow if(epsh.ge.fiemax) epsh=0.9999*fiemax damp=epsh/(fiemax-epsh) damp=dexp(-damp) unod=unod+uparz uinda=uinda+uparz*damp c write(26,868) epshow,uparz,rou(m,5) 141 continue uitota=uitota+uinda write(6,112) k,uinda,unod write(8,112) k,uinda,unod 149 continue 999 continue if(icry.eq.0) flatt=1.0 if(icry.ne.0) flatt=0.5 ectot=ectot*flatt edtot=edtot*flatt ertot=ertot*flatt totota=ectot+uitota+edtot+ertot if(icry.eq.0) write(6,777) ectot,uitota,edtot,ertot,totota if(icry.eq.0) write(8,777) ectot,uitota,edtot,ertot,totota 777 format(/2x,'PIXEL,ec ep ed er et',6f8.1) if(icry.ne.0) write(6,778) ectot,uitota,edtot,ertot,totota if(icry.ne.0) write(8,778) ectot,uitota,edtot,ertot,totota 778 format(/2x,'PIXEL lattice ec ep ed er et',6f8.1) ffatto=ffatto*flatt ffreto=ffreto*flatt fftoto=fftoto*flatt ffqqto=ffqqto*flatt efftot=fftoto+ffqqto if(icry.eq.0) write(6,779) ffatto,ffreto,fftoto,ffqqto,efftot if(icry.eq.0) write(8,779) ffatto,ffreto,fftoto,ffqqto,efftot 779 format(2x,'FF,e6 er e6r eqq etot',5f8.1) if(icry.ne.0) write(6,780) ffatto,ffreto,fftoto,ffqqto,efftot if(icry.ne.0) write(8,780) ffatto,ffreto,fftoto,ffqqto,efftot 780 format(2x,'FF lattice e6 er e6r eqq etot',5f8.1) c branch back to read another box c close electron density files close (unit=11) close (unit=15) read(9,*) ibra if(ibra.eq.0) go to 997 stop end c c this subroutine prepares surrounding molecules for crystal subroutine precry (cutmi,cutma,cuttop,npe) c include all molecules with distance between c centers of mass less than cutoff angstroms implicit real*8 (a-h,o-z) dimension f1(3,3),f2(3,3),t1(3),t2(3),fbuf(3,3),tbuf(3) dimension p(3,3),f2f1(3,3),f2t1(3),u(3),v(3),pmen(3,3),fbpm(3,3) dimension fkf2f1(3,3),fkf2t1(3),fkt2(3),ptk(3),pptk(3,3) common /cry/ fcry(500,3,3),tcry(500,3) c max 500 molecules in crystal cluster npemax=500 c read basic data read(9,*) a,b,c,al,be,ga write(8,25) a,b,c,al,be,ga write(6,25) a,b,c,al,be,ga 25 format(2x,'cell parameters',3f8.3,3f8.2) call matort (a,b,c,al,be,ga,p) call inv(p,pmen) read(9,*) ((f1(i,j),j=1,3),i=1,3) read(9,*) (t1(i),i=1,3) read(9,*) ((f2(i,j),j=1,3),i=1,3) read(9,*) (t2(i),i=1,3) call prma(f2,f1,f2f1) call prtr(f2,t1,f2t1) c read independent space group matrices c identity must be first matrix read(9,*) npz 4 format(10i5) npe=0 do 2 mm=1,npz read(9,*) ((fbuf(j,k),k=1,3),j=1,3) read(9,*) (tbuf(j),j=1,3) c form product P Mk P(-1) call prma(fbuf,pmen,fbpm) call prma(p,fbpm,pptk) c now start coordinate transformations call prma(pptk,f2f1,fkf2f1) call prtr(pptk,f2t1,fkf2t1) call prtr(pptk,t2,fkt2) call prtr(p,tbuf,ptk) if(mm.ne.1.and.cutma.gt.0.01) go to 22 c insert in crystal cluster the input molecules only npe=npe+1 do 23 i=1,3 tcry(npe,i)=fkf2t1(i)+fkt2(i)+ptk(i) do 23 j=1,3 23 fcry(npe,i,j)=fkf2f1(i,j) if(cutma.gt.0.01) go to 22 if(npe.eq.1) go to 2 if (npe.gt.npemax) write(6,99) if(npe.gt.npemax) stop d=(tcry(npe,1)-tcry(1,1))**2 d=d+(tcry(npe,2)-tcry(1,2))**2 d=d+(tcry(npe,3)-tcry(1,3))**2 d=dsqrt(d) write(16,26) npe,d,((fbuf(ii,jj),jj=1,3),ii=1,3), 1 tbuf(1),tbuf(2),tbuf(3) go to 2 c search surrounding molecules within cutoff 22 continue kx=cuttop/a ky=cuttop/b kz=cuttop/c kkx=2*(kx+1)+1 kky=2*(ky+1)+1 kkz=2*(kz+1)+1 do 11 i=1,kkx do 12 j=1,kky do 13 k=1,kkz idx=i-kx-2 idy=j-ky-2 idz=k-kz-2 u(1)=float(idx)+tbuf(1) u(2)=float(idy)+tbuf(2) u(3)=float(idz)+tbuf(3) u1=u(1) u2=u(2) u3=u(3) call prtr(p,u,v) u(1)=fkf2t1(1)+fkt2(1)+v(1) u(2)=fkf2t1(2)+fkt2(2)+v(2) u(3)=fkf2t1(3)+fkt2(3)+v(3) d=(u(1)-tcry(1,1))**2 d=d+(u(2)-tcry(1,2))**2 d=d+(u(3)-tcry(1,3))**2 d=sqrt(d) if(mm.eq.1.and.d.lt.0.01) go to 13 c if(idx.eq.0.and.idy.eq.0.and.idz.eq.0) go to 133 if(d.lt.cutmi.or.d.gt.cutma) go to 13 c molecule is within cutoff or is one of the independent set 133 npe=npe+1 if (npe.gt.npemax) write(6,99) 99 format(2x,'too many molecules, reduce cutoff') if(npe.gt.npemax) stop do 7 ii=1,3 tcry(npe,ii)=u(ii) do 7 jj=1,3 7 fcry(npe,ii,jj)=fkf2f1(ii,jj) write(16,26) npe,d,((fbuf(ii,jj),jj=1,3),ii=1,3), 1 u1,u2,u3 26 format(i6,f7.2,9f5.1,3f8.4) 13 continue 12 continue 11 continue 2 continue write(16,26) 999 return end subroutine prma(a,b,c) real*8 a,b,c dimension a(3,3),b(3,3),c(3,3) do 1 i=1,3 do 1 j=1,3 c(i,j)=0.0 do 3 k=1,3 3 c(i,j)=c(i,j)+a(i,k)*b(k,j) 1 continue return end subroutine prtr(a,b,c) real*8 a,b,c dimension a(3,3),b(3),c(3) do 1 i=1,3 c(i)=0.0 do 3 k=1,3 3 c(i)=c(i)+a(i,k)*b(k) 1 continue return end SUBROUTINE MATORT (A,B,C,ALF,BET,GAM,rm) C COMPUTE ORTHOGONALIZATION MATRIX) implicit real*8 (a-h,o-z) dimension rm(3,3) PZ=0.01745329 CALF=COS(ALF*PZ) SALF=SIN(ALF*PZ) CBET=COS(BET*PZ) SBET=SIN(BET*PZ) CGAM=COS(GAM*PZ) SGAM=SIN(GAM*PZ) RM(1,1)=A*SGAM RM(1,2)=0.0 RM(1,3)=C/SGAM RM(2,1)=A*CGAM RM(2,2)=B RM(2,3)=CALF RM(3,1)=0.0 RM(3,2)=0.0 RM(3,3)=1.-CGAM*CGAM-CBET*CBET-RM(2,3)**2 1 +2.*CGAM*CBET*RM(2,3) RM(3,3)=RM(1,3)*SQRT(RM (3,3)) RM(1,3)=RM(1,3)*(CBET-CGAM*RM(2,3)) RM(2,3)=RM(2,3)*C return end SUBROUTINE euler(flt,flf,flc,eu,pz) implicit real*8 (a-h,o-z) dimension eu(3,3) c1=dcos(flt*pz) s1=dsin(flt*pz) c2=dcos(flf*pz) s2=dsin(flf*pz) c3=dcos(flc*pz) s3=dsin(flc*pz) eu(1,1)=C2*C3 eu(1,2)=-C1*S3+S1*S2*C3 eu(1,3)=S1*S3+C1*S2*C3 eu(2,1)=C2*S3 eu(2,2)=C1*C3+S1*S2*S3 eu(2,3)=-S1*C3+C1*S2*S3 eu(3,1)=-S2 eu(3,2)=S1*C2 eu(3,3)=C1*C2 5 RETURN END SUBROUTINE INV(aa,bb) C INVERT MATRIX aa, STORE RESULT IN bb implicit real*8 (a-h,o-z) dimension aa(3,3),bb(3,3) R11=aa(1,1) R12=aa(1,2) R13=aa(1,3) R21=aa(2,1) R22=aa(2,2) R23=aa(2,3) R31=aa(3,1) R32=aa(3,2) R33=aa(3,3) DET=R11*R22*R33+R12*R23*R31+R21*R32*R13-R13*R22*R31 DET=DET-R11*R23*R32-R12*R21*R33 bb(1,1)=(R22*R33-R23*R32)/DET bb(2,1)=-(R21*R33-R23*R31)/DET bb(3,1)=(R21*R32-R22*R31)/DET bb(1,2)=-(R12*R33-R13*R32)/DET bb(2,2)=(R11*R33-R13*R31)/DET bb(3,2)=-(R11*R32-R12*R31)/DET bb(1,3)=(R12*R23-R13*R22)/DET bb(2,3)=-(R11*R23-R13*R21)/DET bb(3,3)=(R11*R22-R12*R21)/DET RETURN END c c this subroutine reads a gaussian electron density file c and condenses every nrd steps c subroutine sumden(roo,nx,ny,nz,nrd,dx,dy,dz) implicit real*8 (a-h,o-z) dimension rored(150,150),bufde(10,500),rofin(150,150) dimension roo(5000000) rewind 21 rewind 14 nor=500 nre=150 bohr=0.52917 if(nx.le.nor.and.ny.le.nor.and.nz.le.nor) go to 1 write(6,2) nx,ny,nz stop 2 format(2x,'original ro grid too large',3i5) 1 nxs=nx/nrd nys=ny/nrd nzs=nz/nrd if(nxs.lt.nre.and.nys.lt.nre.and.nzs.lt.nre) go to 333 write(6,4) nxs,nys,nzs 4 format(2x,'reduced ro grid too large',3i5) stop 333 nxr=mod(nx,nrd) nyr=mod(ny,nrd) nzr=mod(nz,nrd) c prepare slices of reduced ro density vcel=dx*dy*dz do 101 i=1,nx do 12 j=1,nys do 13 k=1,nrd kre=(j-1)*nrd+k do 137 kk=1,nz iadd=(i-1)*nz*ny+(kre-1)*nz+kk 137 bufde(k,kk)=roo(iadd) c read(ntape) (bufde(k,kk),kk=1,nz) 13 continue 3 format(6e13.5) do 14 l=1,nzs rored(j,l)=0.0 li=(l-1)*nrd+1 lm=l*nrd do 15 kk=1,nrd do 15 m=li,lm 15 rored(j,l)=rored(j,l)+bufde(kk,m)*vcel 14 continue nz1=nzs if(nzr.eq.0) go to 12 c take care of extra points at end of line along z nz1=nzs+1 lm=lm+1 rored(j,nz1)=0.0 do 25 kk=1,nrd do 25 m=lm,nz 25 rored(j,nz1)=rored(j,nz1)+bufde(kk,m)*vcel 12 continue ny1=nys if(nyr.eq.0) go to 37 c take care of extra points at end of line along y ny1=nys+1 do 31 k=1,nyr kre=kre+1 do 6673 kk=1,nz iadd=(i-1)*nz*ny+(kre-1)*nz+kk bufde(k,kk)=roo(iadd) c read(ntape) (bufde(k,kk),kk=1,nz) 6673 continue 31 continue do 34 l=1,nzs li=(l-1)*nrd+1 lm=l*nrd rored(ny1,l)=0.0 do 35 kk=1,nyr do 35 m=li,lm 35 rored(ny1,l)=rored(ny1,l)+bufde(kk,m)*vcel 34 continue if(nyr.eq.0) go to 37 c take care of extra point at end of lines along z and y lm=lm+1 rored(ny1,nz1)=0.0 do 36 kk=1,nyr do 36 m=lm,nz 36 rored(ny1,nz1)=rored(ny1,nz1)+bufde(kk,m)*vcel c end of one layer, write reduced density to disk file 37 continue write(21) ((rored(ik,jk),jk=1,nz1),ik=1,ny1) c close loop on all layers, coordinate along x 101 continue rewind 21 c now start summing up layers along x do 102 i=1,nxs do 111 ii=1,ny1 do 111 jj=1,nz1 111 rofin(ii,jj)=0.0 do 103 j=1,nrd read(21) ((rored(ii,jj),jj=1,nz1),ii=1,ny1) do 104 k=1,ny1 do 104 l=1,nz1 104 rofin(k,l)=rofin(k,l)+rored(k,l) 103 continue do 112 ii=1,ny1 write(14,3) (rofin(ii,jj),jj=1,nz1) 112 continue 102 continue nx1=nxs if(nxr.eq.0) go to 201 c now take care of extra layers at end of x direction nx1=nxs+1 do 107 i=1,ny1 do 107 j=1,nz1 107 rofin(i,j)=0.0 do 105 i=1,nxr read(21) ((rored(ii,jj),jj=1,nz1),ii=1,ny1) do 106 j=1,ny1 do 106 k=1,nz1 106 rofin(j,k)=rofin(j,k)+rored(j,k) 105 continue do 202 ii=1,ny1 write(14,3) (rofin(ii,jj),jj=1,nz1) 202 continue c end business, define new number of points 201 continue write(8,204) nrd, nx,ny,nz,nx1,ny1,nz1 write(6,204) nrd, nx,ny,nz,nx1,ny1,nz1 204 format(' == N.steps,original, condensed',i4,6i5) nx=nx1 ny=ny1 nz=nz1 rewind 14 return end C THIS SUBROUTINE ASSIGNS UNI POTENTIAL ENERGY PARAMETERS subroutine uni6 implicit real*8 (a-h,o-z) character*2 flab COMMON /A/ tipo(30,30,3),rinter(30),rintm(30,30),flab(30) write(8,88) write(6,88) 88 format(' == 6-exp UNI potentials used') c $$$$$$$$$$$$$$$$$$$$$$$ c note: units of kcal/mol c then transformed into kJ/mol c $$$$$$$$$$$$$$$$$$$$$$$ DO 1 I=1,30 DO 1 J=1,30 DO 1 K=1,3 1 TIPO(I,J,K)=0.0 c hydrogen TIPO(1,1,1)=5774.0 TIPO(1,1,2)=4.01 TIPO(1,1,3)=26.1 TIPO(1,3,1)=28870. TIPO(1,3,2)=4.10 TIPO(1,3,3)=113. TIPO(1,4,1)=54560. TIPO(1,4,2)=4.52 TIPO(1,4,3)=120. TIPO(1,5,1)=70610. TIPO(1,5,2)=4.82 TIPO(1,5,3)=105. TIPO(1,6,1)=15358.0 TIPO(1,6,2)=4.12 TIPO(1,6,3)=59.4 TIPO(1,9,1)=64190. TIPO(1,9,2)=4.03 TIPO(1,9,3)=279. TIPO(1,10,1)=70020. TIPO(1,10,2)=4.09 TIPO(1,10,3)=279. TIPO(1,16,1)=70610. TIPO(1,16,2)=4.82 TIPO(1,16,3)=105. TIPO(1,19,1)=54560. TIPO(1,19,2)=4.52 TIPO(1,19,3)=120. TIPO(1,20,1)=54560. TIPO(1,20,2)=4.52 TIPO(1,20,3)=120. TIPO(1,21,1)=64190. TIPO(1,21,2)=4.03 TIPO(1,21,3)=279. TIPO(1,22,1)=70610. TIPO(1,22,2)=4.82 TIPO(1,22,3)=105. TIPO(1,25,1)= 5774. TIPO(1,25,2)=4.01 TIPO(1,25,3)=26.1 c carbon TIPO(3,3,1)=54050. TIPO(3,3,2)=3.47 TIPO(3,3,3)=578. TIPO(3,4,1)=117470. TIPO(3,4,2)=3.86 TIPO(3,4,3)=667. TIPO(3,5,1)=93950. TIPO(3,5,2)=3.74 TIPO(3,5,3)=641. TIPO(3,9,1)=126460. TIPO(3,9,2)=3.41 TIPO(3,9,3)=1504. TIPO(3,10,1)=93370. TIPO(3,10,2)=3.52 TIPO(3,10,3)=923. TIPO(3,16,1)=93950. TIPO(3,16,2)=3.74 TIPO(3,16,3)=641. TIPO(3,19,1)=117470. TIPO(3,19,2)=3.86 TIPO(3,19,3)=667. TIPO(3,20,1)=117470. TIPO(3,20,2)=3.86 TIPO(3,20,3)=667. TIPO(3,21,1)=126460. TIPO(3,21,2)=3.41 TIPO(3,21,3)=1504. TIPO(3,22,1)=93950. TIPO(3,22,2)=3.74 TIPO(3,22,3)=641. TIPO(3,25,1)=28870. TIPO(3,25,2)=4.10 TIPO(3,25,3)=113. c nitrogen TIPO(4,4,1)=87300. TIPO(4,4,2)=3.65 TIPO(4,4,3)=691. TIPO(4,5,1)=64190. TIPO(4,5,2)=3.86 TIPO(4,5,3)=364. TIPO(4,16,1)=64190. TIPO(4,16,2)=3.86 TIPO(4,16,3)=364. TIPO(4,22,1)=64190. TIPO(4,22,2)=3.86 TIPO(4,22,3)=364. TIPO(4,25,1)=5704430. TIPO(4,25,2)=7.78 TIPO(4,25,3)=377. c oxygen TIPO(5,5,1)=46680. TIPO(5,5,2)=3.74 TIPO(5,5,3)=319. TIPO(5,9,1)=110160. TIPO(5,9,2)=3.63 TIPO(5,9,3)=906. TIPO(5,10,1)=80855. TIPO(5,10,2)=3.63 TIPO(5,10,3)=665. TIPO(5,19,1)=64190. TIPO(5,19,2)=3.86 TIPO(5,19,3)=364. TIPO(5,20,1)=64190. TIPO(5,20,2)=3.86 TIPO(5,20,3)=364. TIPO(5,25,1)=70610. TIPO(5,25,2)=4.82 TIPO(5,25,3)=105. c fluorine TIPO(6,6,1)=40850. TIPO(6,6,2)=4.22 TIPO(6,6,3)=135. c silicon TIPO(7,7,1)=232473. TIPO(7,7,2)=3.21 TIPO(7,7,3)=3953. c sulfur TIPO(9,9,1)=259960. TIPO(9,9,2)=3.52 TIPO(9,9,3)=2571. TIPO(9,25,1)=64190. TIPO(9,25,2)=4.03 TIPO(9,25,3)=279. c chlorine TIPO(10,10,1)=140050. TIPO(10,10,2)=3.52 TIPO(10,10,3)=1385. TIPO(10,25,1)=70020. TIPO(10,25,2)=4.09 TIPO(10,25,3)=279. c argon tipo(11,11,1)=164572.0 tipo(11,11,2)=3.553 tipo(11,11,3)=1528.5 c bromine TIPO(13,13,1)=270600. TIPO(13,13,2)=3.28 TIPO(13,13,3)=3628. c iodine TIPO(15,15,1)=372900. TIPO(15,15,2)=3.03 TIPO(15,15,3)=8373. c h-bonding potentials TIPO(16,16,1)=46680. TIPO(16,16,2)=3.74 TIPO(16,16,3)=319. TIPO(16,19,1)=64190. TIPO(16,19,2)=3.86 TIPO(16,19,3)=364. TIPO(16,20,1)=64190. TIPO(16,20,2)=3.86 TIPO(16,20,3)=364. TIPO(16,25,1)=4509750. TIPO(16,25,2)=7.78 TIPO(16,25,3)=298. TIPO(17,17,1)=46680. TIPO(17,17,2)=3.74 TIPO(17,17,3)=319. TIPO(17,19,1)=64190. TIPO(17,19,2)=3.86 TIPO(17,19,3)=364. TIPO(17,20,1)=64190. TIPO(17,20,2)=3.86 TIPO(17,20,3)=364. TIPO(17,25,1)=5336020. TIPO(17,25,2)=8.27 TIPO(17,25,3)=247. TIPO(17,26,1)=6313670. TIPO(17,26,2)=8.75 TIPO(17,26,3)=205. TIPO(17,27,1)=3607810. TIPO(17,27,2)=7.78 TIPO(17,27,3)=238. TIPO(17,28,1)=3607810. TIPO(17,28,2)=7.78 TIPO(17,28,3)=238. TIPO(18,18,1)=46680. TIPO(18,18,2)=3.74 TIPO(18,18,3)=319. TIPO(18,19,1)=64190. TIPO(18,19,2)=3.86 TIPO(18,19,3)=364. TIPO(18,20,1)=64190. TIPO(18,20,2)=3.86 TIPO(18,20,3)=364. TIPO(18,25,1)=4509750. TIPO(18,25,2)=7.78 TIPO(18,25,3)=298. TIPO(18,26,1)=70610. TIPO(18,26,2)=4.82 TIPO(18,26,3)=105. TIPO(18,27,1)=4509750. TIPO(18,27,2)=7.78 TIPO(18,27,3)=298. TIPO(18,28,1)=4509750. TIPO(18,28,2)=7.78 TIPO(18,28,3)=298. TIPO(19,19,1)=87300. TIPO(19,19,2)=3.65 TIPO(19,19,3)=691. TIPO(19,20,1)=87300. TIPO(19,20,2)=3.65 TIPO(19,20,3)=691. TIPO(19,22,1)=64190. TIPO(19,22,2)=3.86 TIPO(19,22,3)=364. TIPO(19,25,1)=5704430. TIPO(19,25,2)=7.78 TIPO(19,25,3)=377. TIPO(19,26,1)=5704430. TIPO(19,26,2)=7.78 TIPO(19,26,3)=377. TIPO(19,27,1)=54560. TIPO(19,27,2)=4.52 TIPO(19,27,3)=120. TIPO(19,28,1)=54560. TIPO(19,28,2)=4.52 TIPO(19,28,3)=120. TIPO(20,20,1)=87300. TIPO(20,20,2)=3.65 TIPO(20,20,3)=691. TIPO(20,22,1)=64190. TIPO(20,22,2)=3.86 TIPO(20,22,3)=364. TIPO(20,25,1)=2852230. TIPO(20,25,2)=7.58 TIPO(20,25,3)=222. TIPO(20,26,1)=2852230. TIPO(20,26,2)=7.58 TIPO(20,26,3)=222. TIPO(20,27,1)=54560. TIPO(20,27,2)=4.52 TIPO(20,27,3)=120. TIPO(20,28,1)=54560. TIPO(20,28,2)=4.52 TIPO(20,28,3)=120. TIPO(21,21,1)=259960. TIPO(21,21,2)=3.52 TIPO(21,21,3)=2571. TIPO(21,25,1)=1750600. TIPO(21,25,2)=5.63 TIPO(21,25,3)=1032. TIPO(22,22,1)=46680. TIPO(22,22,2)=3.74 TIPO(22,22,3)=319. TIPO(22,25,1)=1315000. TIPO(22,25,2)=7.78 TIPO(22,25,3)=86.7 c $$$$$$$$$$$$$$$$$$$$$$$ c note: units of kcal/mol c then transformed into kJ/mol c $$$$$$$$$$$$$$$$$$$$$$$ do 10 i=25,28 do 10 j=1,3 10 tipo(i,i,j)=tipo(1,1,j) do 7 i=1,16 if(i.eq.4) go to 7 do 77 j=26,28 do 77 k=1,3 77 tipo(i,j,k)=tipo(i,25,k) 7 continue do 8 i=1,4 do 8 j=17,18 do 8 k=1,3 8 tipo(i,j,k)=tipo(i,16,k) do 9 i=21,22 do 9 j=26,28 do 9 k=1,3 9 tipo(i,j,k)=tipo(i,25,k) DO 31 Is1=1,30 DO 32 is2=Is1,30 ISTIP=TIPO(Is1,is2,1) IF(ISTIP.GT.0) GO TO 32 TIPO(IS1,IS2,1)=(TIPO(IS1,IS1,1)*TIPO(IS2,IS2,1))**0.5 TIPO(IS1,IS2,3)=(TIPO(IS1,IS1,3)*TIPO(IS2,IS2,3))**0.5 TIPO(IS1,IS2,2)=(TIPO(IS1,IS1,2)+TIPO(IS2,IS2,2))*0.5 32 CONTINUE 31 CONTINUE C SYMMETRIZE POTENTIAL PARAMETER MATRIX DO 2 I=1,30 DO 2 J=I,30 tipo(i,j,1)=tipo(i,j,1)*4.184 tipo(i,j,3)=tipo(i,j,3)*4.184 DO 2 K=1,3 2 TIPO(J,I,K)=TIPO(I,J,K) c load intermolecular radii table do 22 i=1,30 do 22 j=1,30 rintm(i,j)=rinter(i)+rinter(j) 22 continue do 23 i=25,28 rintm(i,4)=2.0 rintm(4,i)=2.0 rintm(i,5)=1.9 rintm(5,i)=1.9 rintm(i,22)=1.9 rintm(22,i)=1.9 rintm(i,21)=1.9 rintm(21,i)=1.9 do 24 j=16,20 rintm(i,j)=1.9 rintm(j,i)=rintm(i,j) 24 continue 23 continue rintm(17,26)=1.7 rintm(26,17)=1.7 RETURN END subroutine doall(iway,nm1,nm2,na1,na2,xa,g1,xb,g2,fnuc1,fnuc2, 1 nr1,ro1,nr2,ro2,iss1,iss2,isp1,isp2,q1,q2,av1,av2,box1,box2, 2 boxmi1,boxmi2,boxma1,boxma2,stp1,stp2,ny1,nz1,ny2,nz2) c calculates electrostatic energy, electrostatic potential, c electrostatic field, exchange repulsion and dispersion, c and point-charge electrostatic energy c also calculates 6-exp intermolecular energies c iway =1 solute-solute c iway =2 solute-solvent c iway =3 solvent-solute c iway =1 solvent-solvent implicit real*8 (a-h,o-z) logical igo, igo3 character*2 flab c f is the electrostatic potential, epse is the electrostatic field dimension epse(60000,3),isp1(50),isp2(50) dimension box1(5000000),box2(5000000),boxmi1(3), 1 boxmi2(3),boxma1(3),boxma2(3) dimension vcub1(8,3),vcub2(8,3),vrot1(8,3),vrot2(8,3) dimension xl1(3),xu1(3),xl2(3),xu2(3),xmint(3),xmant(3) dimension xave1(3),xave2(3) c xa and xb are nuclear coordinates, g1 and g2 are grid molecular c position vectors, fnuc1 and fnuc2 are nuclear charges, c ro1 and ro2 are original ro vectors, rox and roy are rotated c ro vectors, x1 and x2 are rotated nuclear coordinates c c iss1 or iss2 is zero for solute and 1 for solvent molecule c dimension xa(50,3),xb(50,3),g1(500,6),g2(500,6),q1(50),q2(50), 1 fnuc1(50),fnuc2(50),x1(50,3),x2(50,3),ff1(3,3),ff2(3,3),buf(3), 2 ro1(60000,6),ro2(60000,6),rox(60000,3),roy(60000,3) dimension vecin(3),vecro1(3),vecro2(3),buf2(3) dimension ff1m(3,3),ff2m(3,3) COMMON /A/ tipo(30,30,3),rinter(30),rintm(30,30),flab(30) common /b/ eltot,eunir,eunia,eqq,edisp,erep,icry common /ovla/ bdrep,bdexp,collis,ddamp common /cry/ fcry(500,3,3),tcry(500,3) c state constants nro=60000 igo3=iway.eq.3 avoga=6.0221367e+23 bohr=0.52917 coule=1.6021773e-19 pigr=3.1415926536d+00 pz=2*pigr/360.0 epszer=8.85419e-12 angmet=1.0e+10 fourp=4.0*pigr*epszer facpot=coule*angmet/fourp facen=facpot*coule facfie=facpot*angmet eltot=0.0 eunir=0.0 eunia=0.0 eqq=0.0 edisp=0.0 erep=0.0 half=0.5d+00 c start preliminaries for overlap integration stint=0.1*half*(stp1+stp2) stint=stint*stint r1=boxma1(1)-boxmi1(1) r2=boxma1(2)-boxmi1(2) r3=boxma1(3)-boxmi1(3) r4=boxma2(1)-boxmi2(1) r5=boxma2(2)-boxmi2(2) r6=boxma2(3)-boxmi2(3) num1=(r1)/stp1+0.0001 num2=(r2)/stp1+0.0001 num3=(r3)/stp1+0.0001 num4=(r4)/stp2+0.0001 num5=(r5)/stp2+0.0001 num6=(r6)/stp2+0.0001 c write(8,744) num1,num2,num3,num4,num5,num6 744 format(1x,'n.steps in original boxes',3i5,2x,3i5) rmax=dmax1(r1,r2,r3,r4,r5,r6) c 8 corners of first cube vcub1(1,1)=boxmi1(1) vcub1(1,2)=boxmi1(2) vcub1(1,3)=boxmi1(3) vcub1(2,1)=boxmi1(1) vcub1(2,2)=boxma1(2) vcub1(2,3)=boxmi1(3) vcub1(3,1)=boxmi1(1) vcub1(3,2)=boxmi1(2) vcub1(3,3)=boxma1(3) vcub1(4,1)=boxmi1(1) vcub1(4,2)=boxma1(2) vcub1(4,3)=boxma1(3) vcub1(5,1)=boxma1(1) vcub1(5,2)=boxmi1(2) vcub1(5,3)=boxmi1(3) vcub1(6,1)=boxma1(1) vcub1(6,2)=boxma1(2) vcub1(6,3)=boxmi1(3) vcub1(7,1)=boxma1(1) vcub1(7,2)=boxmi1(2) vcub1(7,3)=boxma1(3) vcub1(8,1)=boxma1(1) vcub1(8,2)=boxma1(2) vcub1(8,3)=boxma1(3) c 8 corners of second cube vcub2(1,1)=boxmi2(1) vcub2(1,2)=boxmi2(2) vcub2(1,3)=boxmi2(3) vcub2(2,1)=boxmi2(1) vcub2(2,2)=boxma2(2) vcub2(2,3)=boxmi2(3) vcub2(3,1)=boxmi2(1) vcub2(3,2)=boxmi2(2) vcub2(3,3)=boxma2(3) vcub2(4,1)=boxmi2(1) vcub2(4,2)=boxma2(2) vcub2(4,3)=boxma2(3) vcub2(5,1)=boxma2(1) vcub2(5,2)=boxmi2(2) vcub2(5,3)=boxmi2(3) vcub2(6,1)=boxma2(1) vcub2(6,2)=boxma2(2) vcub2(6,3)=boxmi2(3) vcub2(7,1)=boxma2(1) vcub2(7,2)=boxmi2(2) vcub2(7,3)=boxma2(3) vcub2(8,1)=boxma2(1) vcub2(8,2)=boxma2(2) vcub2(8,3)=boxma2(3) c if not crystal, output header for schakal-type graphics file if(icry.ne.1.and.iway.eq.1) write(27,654) (g1(i,ko),ko=1,6) 654 format('TITL',3f7.3,4f7.1) if(icry.ne.1.and.iway.eq.1) write(27,655) 1.,1.,1.,90.,90.,90. 655 format('CELL',6f8.1) do 311 i=1,nm1 c rotate atoms and ro of i-th molecule c if non-crystal take info from grid points c if crystal take info from space-group matrices c c nuclei if (icry.eq.1) go to 451 call euler (g1(i,4),g1(i,5),g1(i,6),ff1,pz) go to 452 c crystal case 451 do 453 kk=1,3 g1(i,kk)=tcry(i,kk) do 453 ll=1,3 453 ff1(kk,ll)=fcry(i,kk,ll) 452 continue do 1 jj=1,na1 do 2 k=1,3 x1(jj,k)=0.0 do 2 l=1,3 2 x1(jj,k)=x1(jj,k)+ff1(k,l)*xa(jj,l) do 3 l=1,3 3 x1(jj,l)=x1(jj,l)+g1(i,l) ispe=isp1(jj) if(iway.ne.1) go to 1 if(icry.ne.1) write(27,61) flab(ispe),(x1(jj,k),k=1,3) 61 format(2x,'ATOM',2x,a2,3f12.5) 1 continue c ro do 5 jj=1,nr1 do 6 k=1,3 buf(k)=0.0 do 6 l=1,3 6 buf(k)=buf(k)+ff1(k,l)*ro1(jj,l) do 7 l=1,3 7 rox(jj,l)=buf(l)+g1(i,l) 5 continue do 312 j=1,nm2 igo=iway.eq.1.and.i.gt.j if(icry.eq.0) go to 331 c this is a crystalline material c skip all i-j molecule-molecule interactions c not involving molecule number 1 if(i.ne.1.and.j.ne.1) go to 312 331 continue if((iway.eq.1).and.(i.eq.j)) go to 312 c rotate atoms and ro of j-th molecule c nuclei if (icry.eq.1) go to 455 call euler (g2(j,4),g2(j,5),g2(j,6),ff2,pz) go to 456 c crystal case 455 do 457 kk=1,3 g2(j,kk)=tcry(j,kk) do 457 ll=1,3 457 ff2(kk,ll)=fcry(j,kk,ll) 456 continue do 11 jj=1,na2 do 12 k=1,3 x2(jj,k)=0.0 do 12 l=1,3 12 x2(jj,k)=x2(jj,k)+ff2(k,l)*xb(jj,l) do 13 l=1,3 13 x2(jj,l)=x2(jj,l)+g2(j,l) 11 continue c ro do 15 jj=1,nr2 do 16 k=1,3 buf(k)=0.0 do 16 l=1,3 16 buf(k)=buf(k)+ff2(k,l)*ro2(jj,l) do 17 l=1,3 17 roy(jj,l)=buf(l)+g2(j,l) 15 continue c energy, nuclei-nuclei elnn=0.0 if(igo.or.igo3) go to 404 do 403 in=1,na1 do 402 jn=1,na2 rnuc=(x1(in,1)-x2(jn,1))**2+(x1(in,2)-x2(jn,2))**2 rnuc=rnuc+(x1(in,3)-x2(jn,3))**2 rnuc=dsqrt(rnuc) isi=isp1(in) isj=isp2(jn) c nucleus-nucleus electrostatic energy elnn=elnn+fnuc1(in)*fnuc2(jn)/rnuc c point-charge coulomb energy eqq=eqq+q1(in)*q2(jn)/rnuc c check very short internuclear distance rsr=rnuc/(rintm(isi,isj)) if(rsr.gt.0.75) go to 458 write(6,459) in,jn,rnuc write(8,459) in,jn,rnuc 459 format(2x,' ==dist.below 75% Rinter,atoms',2i4,f10.3) if(rnuc.gt.1.0) go to 458 eunir=eunir+10000. eunia=0.0 go to 461 c 6-exp potential 458 eunir=eunir+tipo(isi,isj,1)*exp(-tipo(isi,isj,2)*rnuc) eunia=eunia-tipo(isi,isj,3)*rnuc**(-6) c core repulsion and short-range energy corrections 461 continue 402 continue 403 continue 404 continue do 320 k=1,nr2 do 320 n=1,3 320 epse(k,n)=0.0 elen=0.0 elee=0.0 edi=0.0 edir=0.0 novfi=0 nove=0 volovl=0.0 rosum=0.0 do 313 j1=1,nr2 cha2=ro2(j1,4) xx2=roy(j1,1) yy2=roy(j1,2) zz2=roy(j1,3) c electron-nuclei part do 401 in=1,na1 xn=xx2-x1(in,1) yn=yy2-x1(in,2) zn=zz2-x1(in,3) rnn=xn**2+yn**2+zn**2 rn=dsqrt(rnn) 817 rna=rn if(rn.gt.collis) go to 871 rn=collis novfi=novfi+1 871 rn3=fnuc1(in)/rn**3 c field epse(j1,1)=epse(j1,1)+xn*rn3 epse(j1,2)=epse(j1,2)+yn*rn3 epse(j1,3)=epse(j1,3)+zn*rn3 if(rna.gt.collis) go to 872 rna=collis nove=nove+1 c energy, electrons-nuclei 872 elen=elen+cha2*fnuc1(in)/rna 401 continue polaa=ro2(j1,5) wa=ro2(j1,6) c electron-electron part do 316 i1=1,nr1 cha1=ro1(i1,4) xx=xx2-rox(i1,1) yy=yy2-rox(i1,2) zz=zz2-rox(i1,3) ri=xx**2+yy**2+zz**2 ri=dsqrt(ri) ria=ri rib=ri if(ri.gt.collis) go to 873 ri=collis novfi=novfi+1 873 ri3=cha1/ri**3 epse(j1,1)=epse(j1,1)+xx*ri3 epse(j1,2)=epse(j1,2)+yy*ri3 epse(j1,3)=epse(j1,3)+zz*ri3 if(igo.or.igo3) go to 316 if(ria.gt.collis) go to 874 ria=collis nove=nove+1 874 elee=elee+cha2*cha1/ria 361 continue c dispersion contribution damp=1.0 if(rib.gt.ddamp) go to 460 damp=(ddamp/rib)-1.0 damp=exp(-damp*damp) 460 edi=edi+damp*ro1(i1,5)*polaa*rib**(-6) 316 continue 313 continue if(igo.or.igo3) go to 510 c now rotate extremes of original density boxes c then find extremes of rotated ro's do 585 nn=1,8 do 502 jj=1,3 buf(jj)=0.0 do 502 kk=1,3 502 buf(jj)=buf(jj)+ff1(jj,kk)*vcub1(nn,kk) do 503 kk=1,3 503 vrot1(nn,kk)=buf(kk)+g1(i,kk) 585 continue c now rotate extremes of original density boxes c then find extremes of rotated ro's do 586 nn=1,8 do 504 jj=1,3 buf(jj)=0.0 do 504 kk=1,3 504 buf(jj)=buf(jj)+ff2(jj,kk)*vcub2(nn,kk) do 505 kk=1,3 505 vrot2(nn,kk)=buf(kk)+g2(j,kk) 586 continue c do 1771 kk=1,8 c write(8,1770) (vcub1(kk,kkk),kkk=1,3) c write(8,1772) (vrot1(kk,kkk),kkk=1,3) c 1770 format(10x,3f8.3) c 1772 format(20x,3f8.2) c 1771 continue c do 1773 kk=1,8 c write(8,1770) (vcub2(kk,kkk),kkk=1,3) c write(8,1772) (vrot2(kk,kkk),kkk=1,3) c 1773 continue do 506 kk=1,3 xu1(kk)=-999999. xl1(kk)=999999. xu2(kk)=-999999. 506 xl2(kk)=999999. do 507 nn=1,3 do 507 kk=1,8 if(vrot1(kk,nn).ge.xu1(nn)) xu1(nn)=vrot1(kk,nn) if(vrot2(kk,nn).ge.xu2(nn)) xu2(nn)=vrot2(kk,nn) if(vrot1(kk,nn).le.xl1(nn)) xl1(nn)=vrot1(kk,nn) if(vrot2(kk,nn).le.xl2(nn)) xl2(nn)=vrot2(kk,nn) 507 continue xave1(1)=0.5*(xu1(1)+xl1(1)) xave1(2)=0.5*(xu1(2)+xl1(2)) xave1(3)=0.5*(xu1(3)+xl1(3)) xave2(1)=0.5*(xu2(1)+xl2(1)) xave2(2)=0.5*(xu2(2)+xl2(2)) xave2(3)=0.5*(xu2(3)+xl2(3)) c write(8,511) (xu1(n),xl1(n),xave1(n),n=1,3) c write(8,511) (xu2(n),xl2(n),xave2(n),n=1,3) c integration limits are between higher of minimum coordinates c and lower of maximum coordinates do 509 jj=1,3 xmint(jj)=xl1(jj) if(xl2(jj).gt.xmint(jj)) xmint(jj)=xl2(jj) xmant(jj)=xu1(jj) if(xu2(jj).lt.xmant(jj)) xmant(jj)=xu2(jj) 509 continue c calculate integral overlap c prepare integration limits and step rinte=0.5*(stp1+stp2) nintx=0.00001+(xmant(1)-xmint(1))/rinte ninty=0.00001+(xmant(2)-xmint(2))/rinte nintz=0.00001+(xmant(3)-xmint(3))/rinte rlim=(xave1(1)-xave2(1))**2 rlim=rlim+(xave1(2)-xave2(2))**2 rlim=rlim+(xave1(3)-xave2(3))**2 rlim=dsqrt(rlim) c write(8,512) rlim,rmax c write(8,513) (xmint(kx),xmant(kx),kx=1,3) c write(8,514) nintx,ninty,nintz,rinte 514 format(1x,'n.of integr.pts. and step',3i5,1x,f9.3) 511 format(1x,' rotated box up,low,ave coord'/2x,9f8.3) 512 format(1x,' distance between boxes and max box dimens.'/ 1 1x,2f9.3) 513 format(1x,'integration box dimensions,min-max'/ 1 1x,3(2f9.3,3x)) if(rlim.gt.rmax) go to 510 c prepare vector in integration space vlint=rinte**3 call inv(ff1,ff1m) call inv(ff2,ff2m) do 521 ii=1,nintx vecin(1)=xmint(1)+dfloat(ii-1)*rinte do 521 jj=1,ninty vecin(2)=xmint(2)+dfloat(jj-1)*rinte do 521 kk=1,nintz vecin(3)=xmint(3)+dfloat(kk-1)*rinte vecro1(1)=vecin(1)-g1(i,1) vecro1(2)=vecin(2)-g1(i,2) vecro1(3)=vecin(3)-g1(i,3) vecro2(1)=vecin(1)-g2(j,1) vecro2(2)=vecin(2)-g2(j,2) vecro2(3)=vecin(3)-g2(j,3) do 522 kkk=1,3 buf2(kkk)=0.0 do 522 mmm=1,3 522 buf2(kkk)=buf2(kkk)+ff1m(kkk,mmm)*vecro1(mmm) do 523 kkk=1,3 523 vecro1(kkk)=buf2(kkk) do 524 kkk=1,3 buf2(kkk)=0.0 do 524 mmm=1,3 524 buf2(kkk)=buf2(kkk)+ff2m(kkk,mmm)*vecro2(mmm) do 525 kkk=1,3 525 vecro2(kkk)=buf2(kkk) c now find ro1(vecro1) and ro2(vecro2) do 568 kkk=1,3 if(vecro1(kkk).ge.boxma1(kkk)) go to 521 if(vecro1(kkk).le.boxmi1(kkk)) go to 521 if(vecro2(kkk).ge.boxma2(kkk)) go to 521 if(vecro2(kkk).le.boxmi2(kkk)) go to 521 568 continue i1=1.+(vecro1(1)-boxmi1(1))/stp1 d1=vecro1(1)-(boxmi1(1)+dfloat(i1-1)*stp1) i2=1.+(vecro1(2)-boxmi1(2))/stp1 d2=vecro1(2)-(boxmi1(2)+dfloat(i2-1)*stp1) i3=1.+(vecro1(3)-boxmi1(3))/stp1 d3=vecro1(3)-(boxmi1(3)+dfloat(i3-1)*stp1) i4=1.+(vecro2(1)-boxmi2(1))/stp2 d4=vecro2(1)-(boxmi2(1)+dfloat(i4-1)*stp2) i5=1.+(vecro2(2)-boxmi2(2))/stp2 d5=vecro2(2)-(boxmi2(2)+dfloat(i5-1)*stp2) i6=1.+(vecro2(3)-boxmi2(3))/stp2 d6=vecro2(3)-(boxmi2(3)+dfloat(i6-1)*stp2) c i1,i2,i3 are indexes in box1 c i4,i5,i6 are indexes in box2 i11=i1+1 i21=i2+1 i31=i3+1 d1m=(stp1-d1)**2 d2m=(stp1-d2)**2 d3m=(stp1-d3)**2 d4m=(stp2-d4)**2 d5m=(stp2-d5)**2 d6m=(stp2-d6)**2 d1=d1**2 d2=d2**2 d3=d3**2 d4=d4**2 d5=d5**2 d6=d6**2 iadd=(i1-1)*nz1*ny1+(i2-1)*nz1+i3 r1=box1(iadd) rr1=d1+d2+d3 if(rr1.gt.stint) go to 564 ro11=r1 go to 565 564 continue iadd=(i1-1)*nz1*ny1+(i21-1)*nz1+i3 r2=box1(iadd) rr2=d1+d2m+d3 if(rr2.gt.stint) go to 571 ro11=r2 go to 565 571 continue iadd=(i1-1)*nz1*ny1+(i2-1)*nz1+i31 r3=box1(iadd) rr3=d1+d2+d3m if(rr3.gt.stint) go to 572 ro11=r3 go to 565 572 continue iadd=(i1-1)*nz1*ny1+(i21-1)*nz1+i31 r4=box1(iadd) rr4=d1+d2m+d3m if(rr4.gt.stint) go to 573 ro11=r4 go to 565 573 continue iadd=(i11-1)*nz1*ny1+(i2-1)*nz1+i3 r5=box1(iadd) rr5=d1m+d2+d3 if(rr5.gt.stint) go to 574 ro11=r5 go to 565 574 continue iadd=(i11-1)*nz1*ny1+(i21-1)*nz1+i3 r6=box1(iadd) rr6=d1m+d2m+d3 if(rr6.gt.stint) go to 575 ro11=r6 go to 565 575 continue iadd=(i11-1)*nz1*ny1+(i2-1)*nz1+i31 r7=box1(iadd) rr7=d1m+d2+d3m if(rr7.gt.stint) go to 576 ro11=r7 go to 565 576 continue iadd=(i11-1)*nz1*ny1+(i21-1)*nz1+i31 r8=box1(iadd) rr8=d1m+d2m+d3m if(rr8.gt.stint) go to 577 ro11=r8 go to 565 577 ro11=r1/rr1+r2/rr2+r3/rr3+r4/rr4+r5/rr5+r6/rr6+r7/rr7+r8/rr8 sw1=1./rr1+1./rr2+1./rr3+1./rr4+1./rr5+1./rr6+1./rr7+1./rr8 ro11=ro11/sw1 565 i41=i4+1 i51=i5+1 i61=i6+1 iadd=(i4-1)*ny2*nz2+(i5-1)*nz2+i6 s1=box2(iadd) ss1=d4+d5+d6 if(ss1.gt.stint) go to 566 ro22=s1 go to 567 566 continue iadd=(i4-1)*ny2*nz2+(i51-1)*nz2+i6 s2=box2(iadd) ss2=d4+d5m+d6 if(ss2.gt.stint) go to 591 ro22=s2 go to 567 591 continue iadd=(i4-1)*ny2*nz2+(i5-1)*nz2+i61 s3=box2(iadd) ss3=d4+d5+d6m if(ss3.gt.stint) go to 592 ro22=s3 go to 567 592 continue iadd=(i4-1)*ny2*nz2+(i51-1)*nz2+i61 s4=box2(iadd) ss4=d4+d5m+d6m if(ss4.gt.stint) go to 593 ro22=s4 go to 567 593 continue iadd=(i41-1)*ny2*nz2+(i5-1)*nz2+i6 s5=box2(iadd) ss5=d4m+d5+d6 if(ss5.gt.stint) go to 594 ro22=s5 go to 567 594 continue iadd=(i41-1)*ny2*nz2+(i51-1)*nz2+i6 s6=box2(iadd) ss6=d4m+d5m+d6 if(ss6.gt.stint) go to 595 ro22=s6 go to 567 595 continue iadd=(i41-1)*ny2*nz2+(i5-1)*nz2+i61 s7=box2(iadd) ss7=d4m+d5+d6m if(ss7.gt.stint) go to 596 ro22=s7 go to 567 596 continue iadd=(i41-1)*ny2*nz2+(i51-1)*nz2+i61 s8=box2(iadd) ss8=d4m+d5m+d6m if(ss8.gt.stint) go to 597 ro22=s8 go to 567 597 ro22=s1/ss1+s2/ss2+s3/ss3+s4/ss4+s5/ss5+s6/ss6+s7/ss7+s8/ss8 sw2=1./ss1+1./ss2+1./ss3+1./ss4+1./ss5+1./ss6+1./ss7+1./ss8 ro22=ro22/sw2 567 edir=edir+ro11*ro22*vlint volovl=volovl+vlint rosum=rosum+ro11*ro22 c end of loop over integration steps 521 continue c end of overlap calculation 510 continue c prepare scratch file with i-j molecule electrost. energy,joule c and electrostatic field V/meter elee=elee*facen elen=elen*facen elnn=elnn*facen eltot=eltot+elee+elen+elnn edi=-edi*dsqrt(av1*av2)*1970.0 edisp=edisp+edi edir=bdrep*(edir)**bdexp ea=elee*avoga*0.001 eb=elen*avoga*0.001 ec=elnn*avoga*0.001 write(6,328) i,j,ea,eb,ec,edi,edir,volovl,rosum write(16,329) i,j,ea,eb,ec,edi,edir,novfi,nove, 1 volovl,rosum 328 format(2x,2i4,3f11.1,1x,3f6.1,e11.3) 329 format(2i5,3f12.2,2f7.2,2i5,f6.1,e11.3) erep=erep+edir do 322 j1=1,nr2 do 322 n=1,3 epse(j1,n)=epse(j1,n)*facfie 322 continue c write(12) i,iss1,j,iss2,nr2,elee,elen,elnn write(12) ((epse(k,l),l=1,3),k=1,nr2) c close main loops on molecules 312 continue 311 continue c output stopper of schakal file if(icry.ne.1.and.iway.eq.1) write(27,656) 656 format('END') return end