c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c MIN16 - MINIMISING PROGRAM FOR THE BUCKINGHAM AND FOWLER ELECTROSCTATIC c INTERACTION MODEL TAKEN UP TO THE HEXADECAPOLE LEVEL c The program minimises the energy of interaction between c two distributed multipole arrays carrying c hard repulsive cores by reorienting the carrier molecules. c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c Key features: c c - calculation of electrostatic interaction energy up to quadrupole c - calc of es.int with all terms up to t4 = t.abgd (ie all quadrupole c terms and charge.octopole + dipole.octopole + charge-hexadecapole) c - calc of es.int with all terms up to octopole (t6) c - calc of es.int limited by t6 and hexadecapole c - calc of es.int limited by t7 and hexadecapole c - calc of es.int with all terms up to hexadecapole (t8) c - full minimisation (6 parameters = 3 in translation vector, and c 3 rotation angles of molecule B relative to A) c - limited minimisation - only one of the three translation vector c parameters (R) is minimised c - fixed geometry calculation c - translation vector input in either cartesian or polar form c - input of translation vector from keyboard or file c - output to screen and two types of files c - recovery of molecular multipoles from distributed multipoles c C ver. 25.VII.2004 ----- 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 1990: first version c 10.12.97: creation of output directly compatible with PMIFST (this is the c first modification since 5/93) c 23.11.01: modification to output of molecular multipoles c 25.07.04: index of t_zzzzyyx in ESUM corrected from 64 to 324 by c Anatoliy Volkov, volkov@chem.buffalo.edu c c............................................................................. c c MIN16 can access up to four different files during operation. These are: c c 1. File containing molecular data with user provided name c 2. Optional file containing initial positions of B, user named c 3. Main output file containing results of minimisations, user named c 4. Ancillary output file containing details of individual iterations, c this has a name made up of the characters from the output c file name preceding the dot and extension .MON c c STRUCTURE OF MOLECULAR DATA FILE - This is to contain DMA multipoles for c the two molecules (up to hexadecapole), specification of repulsive radii c on atoms, specification of dummy repulsive centers and various control c parameters: c c LINE 1 - NA number of atoms in A (the fixed molecule) (I4) c LINE 2.. - for i=1,NA x, y, z (in bohr) (free format) c NOTE: dma in spherical tensor notation, in the following order c for each atom (atomic units): 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 LINE .. - repeat the above for NB atoms of B (the displaced c molecule), each x,y,z line can optionally carry an alphanumeric c descriptor of the defined atom, separated from the coordinates c by at least one space c LINE .. - two character, right justified names of atoms in molecules A c and B resp. in the order in which their coordinates and c multipoles were declared. These are for the purpose of c assigning hard sphere radii, names recognised by the program c are given in NAMAT. If a name not in NAMAT is used then the c radius is to be given in the next line c PENULTIMATE LINE - IANGST,IMONIT,NDUMMY,ILEVEL (free format) c IANGST determines whether further input and output will use c distances in Angstr (=1), or bohr (=0) c IMONIT determines whether a brief (=1), intermediate (=2) or c full (=3) .MON listing is required. Brief listing summarises c only the minimised energies, intermediate lists also minimised c geometries, values of selected structural parameters and c attractive and repulsive energies, full c listing gives in addition the progress of c minimisation iteration by iteration and all occurrences of c repulsive surface clashes (in form: numbers of atoms in c an atom pair, sum of their v.d.W. radii,actual separation, and c current value of the repulsive term) c NDUMMY is the number of dummy hard core centres for filling c in parts of repulsive surface not sufficiently defined c by atoms (eg. pi-bonds). If NDUMMY is nonzero then this line c should be followed by NDUMMY lines, each containing x,y,z,r - c the coordinates and radius (in Angstroms) of the additional c centres. The dummy centres are only allowed on molecule A and c up to ten are catered for. c ILEVEL: control parameter for the level of electrostatic c calculation such that c ILEVEL=1, quadrupole level ie. all multipole terms up to c and including quadrupole c ilevel=2, t.abgd level ie. quadrupole level + c.o + c.h + d.o c ilevel=3, octopole level ie. all terms up to and including c octopole c ilevel=4, all terms limited by hexadecapole and t.abgdep ie. c octopole level+c.h+d.h+q.h c ilevel=5, as for ilevel=4 but also the octopole-hexadec term c ilevel=6, hexadecapole level, including the hex-hex term c ilevel=-1, level of calculation to be specified during c program execution c LAST LINE - NSTRP (free format) c where NSTRP defines the number of structural parameters to be c evaluated after minimisation, if it is not zero then NSTRP c lines are to follow each containing: c NPART, NA, NB, NC, ND (free format) c where NPART defines the parameter type and is 1,2,3 for bond c length, bond angle, dihedral angle resp. c NA,NB,NC,ND are atom numbers of which the first 2,3,4 c resp. are required for NPART=1,2,3 but four numbers are c to be put on the line in all cases. A negative number c refers to the relevant dummy repulsive centre c c INPUT OF INITIAL ORIENTATION OF B FROM SCREEN: c c I - selection of computation mode and type of orientation vector: c I=-1 or 1 cartesian translation vector c I=-2 or 2,3 polar translation vector c I<0 calculation for fixed geometry c I>0 minimisation c I=1 or 2 full minimisation c I=3 limited minimisation c P, Q, R /Angstr or bohr, degrees: c if I=+-1 then P,Q,R are cartesian components of c.m. A-B c separation vector, their polar counterparts used for c calculation c if I=+-2 or 3 then P,Q,R are polar coordinates R,Theta,Psi c defining the c.m. separation vector, serving as minimisation C parameters par(1),par(2),par(3) c Theta, Psi, Alpha /degrees: c angles defining orientation of B: Theta,Psi are polar angles c defining orientation of a rotation axis in B.cm cartesian c frame, Alpha is a rotation angle about this axis. These c angles are minimisation parameters par(4),par(5),par(6) resp. c c INPUT OF INITIAL ORIENTATION OF B FROM DISK: c c This should be in the same form as above information with the only c difference that a comment line should precede the two lines of orientation c data. Such three line block ie: c comment line c I, P, Q, R c Theta, Psi, Alpha c can be repated as many times as necessary. The comment line is not used c by the program (can be empty) and is only for legibility of the data file. c If ilevel has been set to -1 then four line blocks are required where c an additional line with the value of ilevel is to follow the comment line. c c............................................................................. c c x0,q0,q1,q2,q3,q4 - input cartesian coordinates and spherical distributed c multipoles for molecules A and B c x0,q0,d0,t0,q0,p0 - initial cartesian coordinates and cartesian distributed c multipole for molecules A and B, these are not altered c in the calculation and are used as reference for all c rotations of B c x,q,d,t,q,p - cartesian coordinates and cartesian distributed c multipoles for A and B following rotation of B c esint - array containing pairwise contributions to c electrostatic energy from various multipole terms, c these are summed into vtot in M/P c reptrm - repulsive contribution from the hard sphere potential c_____________________________________________________________________________ c c implicit real*8 (a-h,o-z) common /blk0/q0(50),q1(50,3),q2(50,5),q3(50,7),q4(50,9), . vdwrad(50),x0(50,3),d0(50,3),t0(50,3,3),u0(50,3,3,3) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit common /hexad/p(50,3,3,3,3),p0(50,3,3,3,3) common /cdum/ndummy,cdummy(10,4) common /mon/imonit,filmon common /blkder/iflag real*8 par(6),r(3),vspar(30),cord(6,3),atmase(19),atmas(50) real*4 vdwaal(19) character filnam*25,coment*76,names(5)*10,spout(30)*22, * filmon*25,filreo*25,LINE*80,LABEL(50)*40 character*2 namat(19),nameat(50),CHRDUM integer*2 nsdef(30,5) data names/'charge','dipole','quadrp.','octop.','hexadec.'/ data namat/' H','He',' C',' N',' O',' F','Ne',' P',' S','Cl', * 'Ar','As','Se','Br','Kr','Sb','Te',' I','Xe'/ data atmase/1.007825d0,4.002603d0,12.d0,14.003074d0,15.994915d0, * 18.998405d0,19.992441d0,30.973765d0,31.972074d0,34.968851d0, * 39.962384d0,74.921596d0,79.916527d0,78.918329d0, * 83.911503d0,120.903816d0,129.906238d0,126.904470d0, * 131.904161d0/ data vdwaal/1.20,1.79,1.7,1.5,1.4,1.35,1.6,1.9,1.85,1.80, * 1.92,2.0,2.0,1.95,1.98,2.2,2.2,2.15,2.18/ c nradii=19 pi=3.1415926535897932d0 conv=pi/180d0 bohr=0.5291771d0 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 c c...Open input and output files, read molecular data c write (*,297) 297 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | MIN16 - MINimisation of up to 16-pole multipole', * ' electrostatic interaction',T79,'|'/ * ' | energy for two molecules', * ' with hard repulsive cores',T79,'|'/ * ' |',76(1H_),'|'/' version 25.VII.2004',T64,'Zbigniew KISIEL'/) c 200 write(*,212)' Name of the data file (with hexadecapoles): ' 212 format(1x/a,$) read(*,'(a)',err=200)filnam open(2,file=filnam,status='old',err=200) c 287 write(*,212)' Name of results file: ' read(*,'(a)',err=287)filnam open(3,file=filnam,status='unknown') c i=0 962 i=i+1 if(filnam(i:i).ne.'.'.and.filnam(i:i).ne.' '.and.filnam(i:i) * .ne.char(0))goto 962 i=i-1 filmon=filnam(1:i) filmon(i+1:i+4)='.mon' write(*, '('' Monitor file: '',a)')filmon open(4,file=filmon,status='unknown') write(3,211)' ' write(4,211)' ' close(4) write(3,199) 199 format(1x,78(1h_)) write(3,290) 290 format(1x/' MINIMISATION OF INTERACTION ENERGY OF MOLECULES A' * ,' AND B, EACH CARRYING'/' ELECTRONIC CHARGE IN DISTRIBUTED' * ,' MULTIPOLE FORM AND A HARD REPULSIVE'/' CORE DEFINED BY ' * ,'Van der Waal''s RADII') write(3,199) c c...number of atoms, comment and multipoles for molecule A c read(2,221)na,coment 221 format(i4,a) write(*,211)coment write(3,211)coment 211 format(1x/a) do 1 i=1,na read(2,*)(x0(i,j),j=1,3) do 225 j=1,3 x(i,j)=x0(i,j) 225 continue C BACKSPACE(2) READ(2,'(A)')LINE NDOTS=0 DO 380 J=1,80 IF(LINE(J:J).EQ.'.')NDOTS=NDOTS+1 IF(NDOTS.EQ.3)GOTO 381 380 CONTINUE 381 DO 382 JJ=J+1,80 IF(LINE(JJ:JJ).EQ.' ')GOTO 383 382 CONTINUE 383 LABEL(I)=LINE(JJ:80) C read(2,*)q0(i) read(2,*)(q1(i,j),j=1,3) read(2,*)(q2(i,j),j=1,5) read(2,*)(q3(i,j),j=1,7) read(2,*)(q4(i,j),j=1,9) write(3,288)i write(3,311)' x,y,z:',(x(i,j),j=1,3),label(i) write(3,289)' charge:',q0(i) write(3,289)' dipole:',(q1(i,j),j=1,3) write(3,289)' qdrple:',(q2(i,j),j=1,5) write(3,310)' octple:',(q3(i,j),j=1,7) write(3,310)' hexdpl:',(q4(i,j),j=1,9) 288 format(' ATOM NUMBER',i3,':') 289 format(a,5f12.6) 311 FORMAT(A,3F12.6,A) 310 format(a,5f12.6/8x,5f12.6) 1 continue c c...number of atoms, comment and multipoles for molecule B c read(2,221)nb,coment write(*,211)coment write(3,211)coment do 2 i=na+1,na+nb read(2,*)(x0(i,j),j=1,3) C BACKSPACE(2) READ(2,'(A)')LINE NDOTS=0 DO 386 J=1,80 IF(LINE(J:J).EQ.'.')NDOTS=NDOTS+1 IF(NDOTS.EQ.3)GOTO 385 386 CONTINUE 385 DO 387 JJ=J+1,80 IF(LINE(JJ:JJ).EQ.' ')GOTO 388 387 CONTINUE 388 LABEL(I)=LINE(JJ:80) C read(2,*)q0(i) read(2,*)(q1(i,j),j=1,3) read(2,*)(q2(i,j),j=1,5) read(2,*)(q3(i,j),j=1,7) read(2,*)(q4(i,j),j=1,9) write(3,288)i write(3,311)' x,y,z:',(x0(i,j),j=1,3),label(i) write(3,289)' charge:',q0(i) write(3,289)' dipole:',(q1(i,j),j=1,3) write(3,289)' qdrple:',(q2(i,j),j=1,5) write(3,310)' octple:',(q3(i,j),j=1,7) write(3,310)' hexdpl:',(q4(i,j),j=1,9) 2 continue c c...names of atoms for the purpose of assigning van der Waal's radii for c the hard sphere repulsion (for A, then for B, in order of cartesians c and multipoles) c write(*,211)' NAMES OF ATOMS AND HARD SPHERE RADII:' write(*,211)' molecule.A' write(3,211)' NAMES OF ATOMS AND HARD SPHERE RADII:' write(3,211)' molecule.A' ramax=0.d0 rbmax=0.d0 do 270 i=1,na+nb if(i.eq.na+1)then write(*,211)' molecule.B' write(3,211)' molecule.B' endif read(2,'(a2)')nameat(i) do 271 j=1,nradii if(nameat(i).eq.namat(j))goto 272 271 continue read(2,*)vdwrad(i) atmas(i)=0.000001d0 CHRDUM=NAMEAT(I) IF(CHRDUM(1:1).EQ.'H')ATMAS(I)=1.007825D0 goto 273 272 vdwrad(i)=vdwaal(j) atmas(i)=atmase(j) 273 write(*,274)nameat(i),vdwrad(i) write(3,274)nameat(i),vdwrad(i) 274 format(10x,a,f15.3) vdwrad(i)=vdwrad(i)/bohr radius=dsqrt(x(i,1)**2+x(i,2)**2+x(i,3)**2)+vdwrad(i) if(i.le.na)then if(radius.gt.ramax)ramax=radius else if(radius.gt.rbmax)rbmax=radius endif 270 continue radius=ramax+rbmax+1.d0 write(3,295)radius*bohr write(*,295)radius*bohr 295 format(1x/' Standoff c.m. distance =',f10.5,' Angstr.') read(2,*)iangst,imonit,ndummy,ilevel if(ilevel.ne.-1.and.(ilevel.lt.1.or.ilevel.gt.6))then write(*,'(1x/'' **** ILEVEL='',i2,'' is out of range''/)') . ilevel stop endif levscr=0 if(ilevel.eq.-1)levscr=1 if(iangst.ne.1.and.iangst.ne.0)iangst=1 if(imonit.lt.1.or.imonit.gt.3)imonit=1 c c...dummy repulsive centres c if(ndummy.eq.0)goto 302 if(ndummy.gt.10)ndummy=10 write(3,211)' DUMMY REPULSIVE CENTRES, coordinates and radii:' do 300 i=1,ndummy read(2,*)(cdummy(i,j),j=1,4) write(3,303)(cdummy(i,j),j=1,4) 303 format(1x,3f12.6,f16.6) do 301 j=1,4 cdummy(i,j)=cdummy(i,j)/bohr 301 continue 300 continue c c...structural parameters to be evaluated c 302 read(2,*,err=1255)nstrp if(nstrp.gt.0.and.nstrp.le.30)then do 1250 j=1,nstrp read(2,*,err=1255)(nsdef(j,i),i=1,5) if(nsdef(j,1).lt.1.or.nsdef(j,1).gt.3)goto 1255 kk=2+nsdef(j,1) do 1256 k=2,kk if(nsdef(j,k).lt.0)then if(-nsdef(j,k).gt.ndummy)goto 1255 goto 1256 endif if(nsdef(j,k).lt.1.or.nsdef(j,k).gt.na+nb)goto 1255 1256 continue if(kk.lt.5)then do 1280 k=kk+1,5 nsdef(j,k)=0 1280 continue endif 1250 continue endif goto 1257 1255 write(*,1258)j,char(7) 1258 format(1x/' ***** ERROR in definition of str. param. number ', * i2,a/) stop c 1257 close(2) c c c...convert spherical to cartesian multipoles c CALL CONVER c c c...optional recovery of molecular multipoles c 500 write(*,212)' Do you want molecular multipoles (1/0) ? ' read(*,501,err=500)i 501 format(i1) if(i.lt.0.or.i.gt.1)goto 500 if(i.eq.1)then CALL RECMUL(1,na,'A') CALL RECMUL(na+1,na+nb,'B') endif close(3) c c c...determine whether initial orientation of B will be input from disk c or from keyboard c nrun=0 251 write(*,212)' Rotation of B from keyboard (=0) or disk (=1) ? ' read(*,'(i5)',err=251)ikbdsk if(ikbdsk.ne.1.and.ikbdsk.ne.0)goto 251 if(ikbdsk.eq.0)goto 252 253 write(*,212)' Name of file containing reorientation data ? ' read(*,'(a)',err=253)filreo open(2,file=filreo,err=253) c c c...beginning of loop for electrostatic energy minimisation for various c starting relative orientations of A and B c c 201 if(ikbdsk.eq.0)goto 252 c c...initial orientation of B from from disk c read(2,'(a)',end=256,err=354)coment 354 if(levscr.eq.1)then read(2,*,end=256,err=256)ilevel if(ilevel.lt.1.or.ilevel.gt.6)goto 256 endif read(2,*,end=256,err=256)iflag,(r(i),i=1,3) if(iflag.eq.555)goto 256 355 read(2,*,err=201)theta,psi,alpha if(iflag.lt.-2.or.iflag.gt.8)goto 201 goto 254 256 close(2) goto 220 c c...initial orientation of B from from keyboard c 252 if(levscr.eq.1)then 950 write(*,'(/'' LEVEL OF CALCULATION (1-6)'', * '', or 555 for EXIT: '',$)') read(*,*,err=950)ilevel if(ilevel.eq.555)goto 220 if(ilevel.lt.1.or.ilevel.gt.6)goto 950 endif write(*,350) 350 format(/' TYPE OF CALCULATION, SPECIFY VALUE FROM -2 to 4:' * /' [-ve - fixed calc only; +ve - fit; 1 - cartes', * 'ian vector AB;', * /' 2,3 - polar vector AB; 3 - translation ', * ' angles theta and psi' * /' fixed and only R,theta1,phi1,alpha fitted;' * /' 4 - fit R only, 555 - exit ]'/' ..... ',$) read(*,*,err=252)iflag if(iflag.eq.555)goto 220 if(iflag.lt.-2.or.iflag.gt.4)goto 252 351 write(*,211)' Vector from A to B (x,y,z or R,theta,phi):' read(*,*,err=351)(r(i),i=1,3) 352 write(*,211)' Rotation of B (angles theta1, phi1, alpha) ? ' read(*,*,err=352)theta,psi,alpha c c...conversion to bohr and of B vector to polar coordinates (if necessaary) c 254 ifit=1 if(iflag.lt.0)then ifit=0 iter=0 iflag=-iflag endif if(iflag.gt.1)then par(1)=r(1) par(2)=r(2)*conv par(3)=r(3)*conv r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) endif c if(iangst.eq.0)goto 210 do 215 i=1,3 r(i)=r(i)/bohr 215 continue c 210 par(1)=dsqrt(r(1)**2+r(2)**2+r(3)**2) if(par(1).eq.0.d0)par(1)=1.d-10 par(2)=r(3)/par(1) if(dabs(par(2)).gt.1.d0)par(2)=dsign(1.d0,par(2)) par(2)=dacos(par(2)) if(r(1).eq.0.0d0)r(1)=1.d-10 par(3)=datan2(r(2),r(1)) par(4)=theta*conv par(5)=psi*conv par(6)=alpha*conv c write(*,199) write(*,281)' ',' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) open(3,file=filnam,access='APPEND') write(3,199) write(3,281)' ',' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) close(3) if(imonit.eq.3)then open(4,file=filmon,access='APPEND') write(4,199) write(4,281)' **********************************************', * ' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) close(4) endif c c c...MINIMISATION c c 291 nrun=nrun+1 if(ifit.eq.1)then CALL DFPMIN(par,6,1.d-8,iter,fret) endif vtot=ESINTE(par) c c c...final results; energy c open(3,file=filnam,access='APPEND') write(3,280)' MINIMISED ENERGY =',vtot,vtot*219476.d0,iter,nrun, * ilevel 280 format(1x/a,f12.8,' <<',f11.2,'cm-1 NIT=',i2,' RUN=',i3, * ' ILEVEL=',i1,'>>') r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) c c...coordinates c write(3,2211)NA+NB, * ' CARTESIAN COORDS. + MASSES FOR THE A-B SYSTEM:' 2211 FORMAT(i5,A) write(*,211)' CARTESIAN COORDINATES FOR THE A-B SYSTEM:' if(iangst.eq.0)goto 117 do 116 i=1,3 r(i)=r(i)*bohr do 116 j=1,na+nb x(j,i)=x(j,i)*bohr 116 continue 117 do 30 i=1,na write(*,206)nameat(i),(x(i,j),j=1,3) write(3,1206)i,(x(i,j),j=1,3),atmas(i),nameat(i) 30 continue do 3099 i=na+1,na+nb write(*,206)nameat(i),x(i,1)+r(1),x(i,2)+r(2),x(i,3)+r(3) write(3,1206)i,x(i,1)+r(1),x(i,2)+r(2),x(i,3)+r(3), * atmas(i),nameat(i) 3099 continue 206 format(1x,a,3f14.6) 1206 format(1x,i2,4f14.6,3x,a) c c...structural parameters c if(nstrp.gt.0)then do 1252 i=1,nstrp do 1253 j=2,5 jj=j-1 if(nsdef(i,1).eq.1.and.j.gt.3)goto 1253 if(nsdef(i,1).eq.2.and.j.gt.4)goto 1253 nn=nsdef(i,j) if(nn.lt.0)then nn=-nn cord(jj,1)=cdummy(nn,1)*bohr cord(jj,2)=cdummy(nn,2)*bohr cord(jj,3)=cdummy(nn,3)*bohr goto 1253 endif if(nn.le.na)then cord(jj,1)=x(nn,1) cord(jj,2)=x(nn,2) cord(jj,3)=x(nn,3) endif if(nn.gt.na)then cord(jj,1)=x(nn,1)+r(1) cord(jj,2)=x(nn,2)+r(2) cord(jj,3)=x(nn,3)+r(3) endif 1253 continue if(nsdef(i,1).eq.1)then vspar(i)=bondl(cord,1,2) write(spout(i),1270)nsdef(i,2),nsdef(i,3),vspar(i) endif if(nsdef(i,1).eq.2)then vspar(i)=angle(cord,conv,1,2,3) write(spout(i),1271)nsdef(i,2),nsdef(i,3), * nsdef(i,4),vspar(i) endif if(nsdef(i,1).eq.3)then angl1=angle(cord,conv,1,2,3)-90.d0 angl2=angle(cord,conv,2,3,4)-90.d0 r12s=bondl(cord,1,2)*dsin(angl1*conv) r23=bondl(cord,2,3) r34s=bondl(cord,3,4)*dsin(angl2*conv) do 1265 k=1,3 u23k=(cord(3,k)-cord(2,k))/r23 cord(5,k)=cord(1,k)-cord(2,k)+r12s*u23k cord(6,k)=cord(4,k)-cord(3,k)-r34s*u23k 1265 continue b1=dsqrt(cord(5,1)**2+cord(5,2)**2+cord(5,3)**2) b2=dsqrt(cord(6,1)**2+cord(6,2)**2+cord(6,3)**2) dotp=cord(5,1)*cord(6,1)+cord(5,2)*cord(6,2) * +cord(5,3)*cord(6,3) vspar(i)=dotp/(b1*b2) if(dabs(vspar(i)).gt.1.d0)vspar(i)=dsign(1.d0,vspar(i)) vspar(i)=dacos(vspar(i))/conv write(spout(i),1272)nsdef(i,2),nsdef(i,3), * nsdef(i,4),nsdef(i,5),vspar(i) endif c 1252 continue c write(*,1273)(spout(i),i=1,nstrp) write(3,1273)(spout(i),i=1,nstrp) if(imonit.gt.1)then open(4,file=filmon,access='APPEND') write(4,1273)(spout(i),i=1,nstrp) close(4) endif 1273 format(1x/10(1x,a22,5x,a22,5x,a22/)) c endif 1270 format('R(',i2,',',i2,')=',f8.5,' ') 1271 format('A(',i2,',',i2,',',i2,')=',f7.2,' ') 1272 format('D(',i2,',',i2,',',i2,',',i2,')=',f7.2) c c...translation vector c if(ifit.eq.1)then write(3,281)' ',' MINIMISED',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) write(*,281)' ',' MINIMISED',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) if(imonit.gt.1)then open(4,file=filmon,access='APPEND') write(4,281)' ----------------------------------------', * ' MINIMISED',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv,(par(i)/conv,i=4,6) close(4) endif endif 281 format(1x,a/a,' DISPLACEMENT VECTOR AND ROTATION OF B:'/ * ' cartesian: ',3f15.6/' polar: ',3f15.6 * /' Theta, Psi, Alpha: ',3f15.6) c if(iangst.eq.0)goto 131 do 132 i=1,3 r(i)=r(i)/bohr do 132 j=1,na+nb x(j,i)=x(j,i)/bohr 132 continue c c...contributions to energy c 131 write(*,282)(names(i),i=1,4) write(3,282)(names(i),i=1,4) 282 format(/' CONTRIBUTIONS TO AND TOTAL INTERACTION ENERGY:'/ * 14x,4(2x,a10),7x,'Total') vtot=0.d0 do 244 i=1,5 rtot=0.d0 do 31 j=1,i vtot=vtot+esint(i,j) rtot=rtot+esint(i,j) 31 continue if(i.eq.5)rtot=rtot-esint(5,5) write(*,245)names(i),(esint(i,j),j=1,4),rtot write(3,245)names(i),(esint(i,j),j=1,4),rtot 244 continue 245 format(1x,a10,4f12.7,f16.7) if(esint(5,5).ne.0.d0)then write(*,400)esint(5,5) write(3,400)esint(5,5) endif 400 format(' hexadec-hexadec term ',36x,f16.7) c1tot=esint(1,1)+esint(2,1)+esint(3,1)+esint(4,1)+esint(5,1) c2tot=esint(2,2)+esint(3,2)+esint(4,2)+esint(5,2) c3tot=esint(3,3)+esint(4,3)+esint(5,3) c4tot=esint(4,4)+esint(5,4) write(*,276)c1tot,c2tot,c3tot,c4tot,vtot write(3,276)c1tot,c2tot,c3tot,c4tot,vtot 276 format(1x/1x,'Total ',4f12.7,f16.7) write(*,275)reptrm,vtot+reptrm write(3,275)reptrm,vtot+reptrm 275 format(12x,' + repulsive contribution of ',f16.7,' =', * f14.7) close(3) c open(4,file=filmon,access='APPEND') write(4,292)' MINIMISED ENERGY =',vtot+reptrm, * (vtot+reptrm)*219476.d0,iter,nrun,ilevel 292 format(a,f12.8,' <<',f11.2,'cm-1 NIT=',i3,' RUN=',i3,', ILEVEL=', * i1,'>>') if(imonit.gt.1)then write(4,591)vtot,reptrm endif close(4) 591 format(' attractive energy =',f19.10, * ' repulsive energy =',f19.10) c goto 201 c c 220 open(3,file=filnam,access='APPEND') write(3,199) open(4,file=filmon,access='APPEND') write(4,199) close(3) close(4) c stop end c c_____________________________________________________________________________ c function angle(c,conv,na,nb,nc) c real*8 c(6,3),conv,b1,b2,b3,angle,bondl c b1=bondl(c,na,nb) b2=bondl(c,nb,nc) b3=bondl(c,na,nc) angle=(b1**2+b2**2-b3**2)/(2.d0*b1*b2) if(dabs(angle).gt.1.d0)angle=dsign(1.d0,angle) angle=dacos(angle)/conv c return end c c_____________________________________________________________________________ c function bondl(c,na,nb) c real*8 c(6,3),bondl c bondl=dsqrt((c(na,1)-c(nb,1))**2+ * (c(na,2)-c(nb,2))**2+(c(na,3)-c(nb,3))**2) c return end c c_____________________________________________________________________________ c SUBROUTINE CONVER c c Conversion from spherical to cartesian multipoles c implicit real*8 (a-h,o-z) common /blk0/q0(50),q1(50,3),q2(50,5),q3(50,7),q4(50,9), . vdwrad(50),x0(50,3),d0(50,3),t0(50,3,3),u0(50,3,3,3) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit common /hexad/p(50,3,3,3,3),p0(50,3,3,3,3) real*8 work(81) integer*2 i,j,k,l,m c st=0.5d0*dsqrt(3.0d0) half=0.5d0 sqr58=dsqrt(5.d0/8.d0) sqr38=dsqrt(3.d0/8.d0) sqr23=dsqrt(2.d0/3.d0) sqr512=dsqrt(5.d0/12.d0) sqr52=dsqrt(5.d0/2.d0) sqr5=dsqrt(5.d0) sqr35=dsqrt(35.d0) sqr352=dsqrt(35.d0/2.d0) one8=0.125d0 one3=1.d0/3.d0 c c do 3 i=1,na+nb c c---charge q(i)=q0(i) c c---dipole d0(i,1)=q1(i,2) d0(i,2)=q1(i,3) d0(i,3)=q1(i,1) c c---quadrupole t0(i,1,1)=st*q2(i,4)-half*q2(i,1) t0(i,1,2)=st*q2(i,5) t0(i,2,1)=st*q2(i,5) t0(i,1,3)=st*q2(i,2) t0(i,3,1)=st*q2(i,2) t0(i,2,2)=-st*q2(i,4)-half*q2(i,1) t0(i,2,3)=st*q2(i,3) t0(i,3,2)=st*q2(i,3) t0(i,3,3)=q2(i,1) c c---octopole work( 1)= sqr58* q3(i,6)- sqr38*q3(i,2) work( 4)=-sqr58* q3(i,6)-one3*sqr38*q3(i,2) work( 9)= sqr23* q3(i,2) work( 6)= sqr512*q3(i,5) work( 8)=-sqr58* q3(i,7)- sqr38*q3(i,3) work( 2)= sqr58* q3(i,7)-one3*sqr38*q3(i,3) work(18)= sqr23* q3(i,3) work(27)= q3(i,1) work( 3)= sqr512*q3(i,4)- 0.5d0*q3(i,1) work(12)=-sqr512*q3(i,4)- 0.5d0*q3(i,1) do 1 ia=1,3 do 1 ib=1,3 do 1 ig=1,3 index=ia*ib*ig u0(i,ia,ib,ig)=work(index) 1 continue c c---hexadecapole work( 1)=one8*(3.d0*q4(i,1)-2.d0*sqr5*q4(i,4)+sqr35*q4(i,8)) pxxxx work( 2)=one8*(sqr35*q4(i,9)-sqr5*q4(i,5)) pxxxy work( 3)=one8*(sqr352*q4(i,6)-3.d0*sqr52*q4(i,2)) pxxxz work( 4)=one8*(q4(i,1)-sqr35*q4(i,8)) pxxyy work( 6)=one8*(sqr352*q4(i,7)-sqr52*q4(i,3)) pxxyz work( 9)=0.25d0*(sqr5*q4(i,4)-2.d0*q4(i,1)) pxxzz work( 8)=-one8*(sqr35*q4(i,9)+sqr5*q4(i,5)) pxyyy work(12)=-one8*(sqr352*q4(i,6)+sqr52*q4(i,2)) pxyyz work(18)=0.25d0*sqr5*q4(i,5) pxyzz work(27)=0.5d0*sqr52*q4(i,2) pxzzz work(16)=one8*(3.d0*q4(i,1)+2.d0*sqr5*q4(i,4)+sqr35*q4(i,8)) pyyyy work(24)=-one8*(sqr352*q4(i,7)+3.d0*sqr52*q4(i,3)) pyyyz work(36)=-0.25d0*(sqr5*q4(i,4)+2.d0*q4(i,1)) pyyzz work(54)=0.5d0*sqr52*q4(i,3) pyzzz work(81)=q4(i,1) pzzzz do 2 ia=1,3 do 2 ib=1,3 do 2 ig=1,3 do 2 id=1,3 index=ia*ib*ig*id p0(i,ia,ib,ig,id)=work(index) 2 continue c c 3 continue c c Transfer the multipoles for A to work arrays d,t,u,p, this is really for c neatness since these are not changed during the calculation c do 230 i=1,na do 230 j=1,3 d(i,j)=d0(i,j) do 230 k=1,3 t(i,j,k)=t0(i,j,k) do 230 l=1,3 u(i,j,k,l)=u0(i,j,k,l) do 230 m=1,3 p(i,j,k,l,m)=p0(i,j,k,l,m) 230 continue c return end c c_____________________________________________________________________________ c SUBROUTINE RECMUL(istart,iend,molnam) c c Calculation of molecular multipole for sites from ISTART to IEND c (eq.11 in A.J.Stone, CPL 83,233-239,1981 provides a general expression c for spherical multipoles up to any order) c C Cartesian molecular multipoles are summed into: C charge - qtot C dipole - dtot(3) C quadrupole - ttot(3,3) C octopole - utot(3,3,3) C implicit real*8 (a-h,o-z) parameter (dconv=2.541766d0, tconv=1.345044d0, uconv=0.7117668d0) common /blk0/q0(50),q1(50,3),q2(50,5),q3(50,7),q4(50,9), . vdwrad(50),x0(50,3),d0(50,3),t0(50,3,3),u0(50,3,3,3) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit character molnam character*3 dipole(3),quadrp(5),octop(7) real*8 out(15) data dipole/'10 ','11c','11s'/ data quadrp/'20 ','21c','21s','22c','22s'/ data octop/'30 ','31c','31s','32c','32s','33c','33s'/ c c...calculation of the Cartesian molecular multipole (ZK notebook R3/112) c one3=1.d0/3.d0 qtot=0.0d0 do 8 i=1,3 dtot(i)=0.0d0 do 8 j=1,3 ttot(i,j)=0.0d0 do 8 k=1,3 utot(i,j,k)=0.d0 8 continue 7 continue c do 9 i=istart,iend asq=x0(i,1)**2+x0(i,2)**2+x0(i,3)**2 qtot=qtot+q(i) c do 10 ia=1,3 dtot(ia)=d0(i,ia)+q(i)*x0(i,ia)+dtot(ia) c do 11 ib=1,3 ttot(ia,ib)=t0(i,ia,ib)+1.5d0*(d0(i,ia)*x0(i,ib) . +d0(i,ib)*x0(i,ia)) . -del(ia,ib)*(d0(i,1)*x0(i,1)+d0(i,2)*x0(i,2) . +d0(i,3)*x0(i,3))+1.5d0*q(i)*x0(i,ia)*x0(i,ib) . -0.5d0*del(ia,ib)*q(i)*asq+ttot(ia,ib) c do 12 ig=1,3 sum1=0.d0 sum2=0.d0 do 13 id=1,3 sum1=sum1+x0(i,id)*(t0(i,ia,id)*del(ib,ig) . +t0(i,ib,id)*del(ia,ig)+t0(i,ig,id)*del(ia,ib)) sum2=sum2+x0(i,id)*d0(i,id) 13 continue term=x0(i,ia)*del(ib,ig) . +x0(i,ib)*del(ia,ig)+x0(i,ig)*del(ia,ib) term1=-one3*(2.d0*sum1-5.d0*(x0(i,ia)*t0(i,ib,ig) . +x0(i,ib)*t0(i,ia,ig)+x0(i,ig)*t0(i,ia,ib))) term2=5.d0*(d0(i,ia)*x0(i,ib)*x0(i,ig) . +d0(i,ib)*x0(i,ia)*x0(i,ig) . +d0(i,ig)*x0(i,ia)*x0(i,ib)) term2=0.5d0*(term2-asq*(d0(i,ia)*del(ib,ig) . +d0(i,ib)*del(ia,ig)+d0(i,ig)*del(ia,ib)) . -2.d0*sum2*term) term3=-0.5d0*q(i)*(asq*term . -5.0d0*x0(i,ia)*x0(i,ib)*x0(i,ig)) utot(ia,ib,ig)=utot(ia,ib,ig)+u0(i,ia,ib,ig) . +term1+term2+term3 12 continue 11 continue 10 continue 9 continue C C...Cartesian molecular multipole output - disk and scree C (in experimental units of the GAMESS molecular C multipole listng) C write(*,202)molnam write(3,202)molnam 202 format(1x/' MOLECULAR MULTIPOLES - Molecule.',a1, . /1x,34(1h-)//' Cartesian:'/) c write(*,210)qtot write(3,210)qtot 210 format(1x, 'Charge: ',f20.6) c write(*,205)'Dipole /D:' write(3,205)'Dipole /D:' write(*,206)'X ',dtot(1)*dconv,'Y ',dtot(2)*dconv, . 'Z ',dtot(3)*dconv write(3,206)'X ',dtot(1)*dconv,'Y ',dtot(2)*dconv, . 'Z ',dtot(3)*dconv c write(*,205)'Quadrupole /DA:' write(3,205)'Quadrupole /DA:' write(*,206)'XX ',ttot(1,1)*tconv,'YY ',ttot(2,2)*tconv, . 'ZZ ',ttot(3,3)*tconv write(*,206)'XY ',ttot(1,2)*tconv,'XZ ',ttot(1,3)*tconv, . 'YZ ',ttot(2,3)*tconv write(3,206)'XX ',ttot(1,1)*tconv,'YY ',ttot(2,2)*tconv, . 'ZZ ',ttot(3,3)*tconv write(3,206)'XY ',ttot(1,2)*tconv,'XZ ',ttot(1,3)*tconv, . 'YZ ',ttot(2,3)*tconv c write(*,205)'Octopole /10**(-34) esu:' write(3,205)'Octopole /10**(-34) esu:' write(*,206)'XXX',utot(1,1,1)*uconv,'XXY',utot(1,1,2)*uconv, . 'XXZ',utot(1,1,3)*uconv write(*,206)'XYY',utot(1,2,2)*uconv,'YYY',utot(2,2,2)*uconv, . 'YYZ',utot(2,2,3)*uconv write(*,206)'XZZ',utot(1,3,3)*uconv,'YZZ',utot(2,3,3)*uconv, . 'ZZZ',utot(3,3,3)*uconv write(*,206)'XYZ',utot(1,2,3)*uconv write(3,206)'XXX',utot(1,1,1)*uconv,'XXY',utot(1,1,2)*uconv, . 'XXZ',utot(1,1,3)*uconv write(3,206)'XYY',utot(1,2,2)*uconv,'YYY',utot(2,2,2)*uconv, . 'YYZ',utot(2,2,3)*uconv write(3,206)'XZZ',utot(1,3,3)*uconv,'YZZ',utot(2,3,3)*uconv, . 'ZZZ',utot(3,3,3)*uconv write(3,206)'XYZ',utot(1,2,3)*uconv C C...spherical molecular multipole and output - disk only C (in atomic units) C write(3,211) 211 format(1x/' Spherical:'/) c c write(*,210)qtot write(3,210)qtot c c write(*,205)'Dipole:' write(3,205)'Dipole:' c write(*,206)dipole(1),dtot(3),dipole(2),dtot(1), c . dipole(3),dtot(2) write(3,206)dipole(1),dtot(3),dipole(2),dtot(1), . dipole(3),dtot(2) c c write(*,205)'Quadrupole:' write(3,205)'Quadrupole:' out(1)=ttot(3,3) out(2)=dsqrt(4.d0/3.d0)*ttot(1,3) out(3)=dsqrt(4.d0/3.d0)*ttot(2,3) out(4)=dsqrt(1.d0/3.d0)*(ttot(1,1)-ttot(2,2)) out(5)=dsqrt(4.d0/3.d0)*ttot(1,2) c write(*,206)(quadrp(i),out(i),i=1,5) write(3,206)(quadrp(i),out(i),i=1,5) c c write(*,205)'Octopole:' write(3,205)'Octopole:' 205 format(1x,a) 206 format(3(1x,a3,f11.6,3(5x,a3,f11.6)/)) out(1)=utot(3,3,3) out(2)=dsqrt(3.d0/2.d0)*utot(1,3,3) out(3)=dsqrt(3.d0/2.d0)*utot(2,3,3) out(4)=dsqrt(3.d0/5.d0)*(utot(1,1,3)-utot(2,2,3)) out(5)=2.d0*dsqrt(3.d0/5.d0)*utot(1,2,3) out(6)=dsqrt(1.d0/10.d0)*(utot(1,1,1)-3.d0*utot(1,2,2)) out(7)=dsqrt(1.d0/10.d0)*(3.d0*utot(1,1,2)-utot(2,2,2)) c write(*,206)(octop(i),out(i),i=1,7) write(3,206)(octop(i),out(i),i=1,7) c return end c c______________________________________________________________________ c FUNCTION ESINTE(par) c c Calculation of interaction energy between molecules A and B being the c sum of electrostatic i.e. of multipole arrays at A and B and a repulsive c term from the hard sphere potential c implicit real*8 (a-h,o-z) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit common /hexad/p(50,3,3,3,3),p0(50,3,3,3,3) real*8 par(6) c c...reorient B c CALL ROTATB(par) c c...evaluate the electrostatic interaction energy c CALL ESUM(par) esinte=esint(1,1)+esint(2,1)+esint(2,2) * +esint(3,1)+esint(3,2)+esint(3,3) * +esint(4,1)+esint(4,2)+esint(4,3)+esint(4,4) * +esint(5,1)+esint(5,2)+esint(5,3)+esint(5,4) * +esint(5,5)+reptrm c return end c c_____________________________________________________________________________ c SUBROUTINE ESINTD(par,grad,f) c c Evaluation of derivatives GRAD of interaction energy w.r.t. the six c minimisation parameters in PAR. F is the interaction energy at c PAR. c implicit real*8 (a-h,o-z) common /blkder/iflag real*8 par(6),grad(6),delta(6) data delta/0.0001d0,5*0.00017d0/ c do 1 n=1,6 if(iflag.eq.4.and.n.gt.1)then grad(n)=0.d0 goto 1 endif par(n)=par(n)+delta(n) fd=esinte(par) grad(n)=(fd-f)/delta(n) par(n)=par(n)-delta(n) 1 continue if(iflag.eq.3)then grad(2)=0.d0 grad(3)=0.d0 endif c return end c c_____________________________________________________________________________ c SUBROUTINE ROTATB(par) c c Transformation of coordinates and multipoles of B using the two polar c angles and the rotation angle, the conical transformation matrix c c ( x = c(xnew,xold) * x0 ) is from "Rotations, Quaternions and double c groups", S.L.Altman, Clarendon Press (1986), p.75 c c PAR - parameters being fitted, par(4),par(5),par(6) are theta,psi,alpha c in radians c implicit real*8 (a-h,o-z) common /blk0/q0(50),q1(50,3),q2(50,5),q3(50,7),q4(50,9), . vdwrad(50),x0(50,3),d0(50,3),t0(50,3,3),u0(50,3,3,3) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit common /hexad/p(50,3,3,3,3),p0(50,3,3,3,3) real*8 par(6) integer*2 i,k,l,m,n,ia,ib,ic,id c CALL NORMAL(1,par(4)) CALL NORMAL(2,par(5)) CALL NORMAL(2,par(6)) theta=par(4) psi=par(5) alpha=par(6) c an(1)=dsin(theta)*dcos(psi) an(2)=dsin(theta)*dsin(psi) an(3)=dcos(theta) s=dsin(alpha) ss=(dsin(0.5d0*alpha))**2 c(1,1)=1.d0 - 2.0d0*(an(2)**2+an(3)**2)*ss c(1,2)=-an(3)*s+2.0d0*an(1)*an(2)*ss c(1,3)= an(2)*s+2.0d0*an(1)*an(3)*ss c(2,1)= an(3)*s+2.0d0*an(1)*an(2)*ss c(2,2)=1.d0 - 2.0d0*(an(1)**2+an(3)**2)*ss c(2,3)=-an(1)*s+2.0d0*an(2)*an(3)*ss c(3,1)=-an(2)*s+2.0d0*an(3)*an(1)*ss c(3,2)= an(1)*s+2.0d0*an(2)*an(3)*ss c(3,3)=1.d0 - 2.0d0*(an(1)**2+an(2)**2)*ss c c c...rotate coordinates and moments of B sites c do 21 i=na+1,na+nb c c---coordinates do 22 k=1,3 x(i,k)=0.0d0 do 22 l=1,3 x(i,k)=c(k,l)*x0(i,l)+x(i,k) 22 continue c c---dipoles do 24 k=1,3 d(i,k)=0.0d0 do 24 l=1,3 d(i,k)=c(k,l)*d0(i,l)+d(i,k) 24 continue c c---quadrupoles do 26 k=1,3 do 26 l=1,3 t(i,k,l)=0.0d0 do 26 m=1,3 do 26 n=1,3 t(i,k,l)=c(k,m)*c(l,n)*t0(i,m,n)+t(i,k,l) 26 continue c c---octopoles do 50 m=1,3 do 50 l=1,m do 50 k=1,l zz=0.d0 do 52 ia=1,3 do 52 ib=1,3 do 52 ic=1,3 zz=c(k,ia)*c(l,ib)*c(m,ic) . *u0(i,ia,ib,ic)+zz 52 continue u(i,k,l,m)=zz u(i,k,m,l)=zz u(i,l,k,m)=zz u(i,l,m,k)=zz u(i,m,k,l)=zz u(i,m,l,k)=zz 50 continue c c---hexadecapoles do 51 n=1,3 do 51 m=1,n do 51 l=1,m do 51 k=1,l zz=0.d0 do 53 ia=1,3 do 53 ib=1,3 do 53 ic=1,3 do 53 id=1,3 zz=c(k,ia)*c(l,ib)*c(m,ic)*c(n,id) . *p0(i,ia,ib,ic,id)+zz 53 continue p(i,k,l,m,n)=zz p(i,k,l,n,m)=zz p(i,k,m,l,n)=zz p(i,k,m,n,l)=zz p(i,k,n,l,m)=zz p(i,k,n,m,l)=zz p(i,l,k,m,n)=zz p(i,l,k,n,m)=zz p(i,l,m,k,n)=zz p(i,l,m,n,k)=zz p(i,l,n,m,k)=zz p(i,l,n,k,m)=zz p(i,m,k,l,n)=zz p(i,m,k,n,l)=zz p(i,m,l,k,n)=zz p(i,m,l,n,k)=zz p(i,m,n,k,l)=zz p(i,m,n,l,k)=zz p(i,n,k,l,m)=zz p(i,n,k,m,l)=zz p(i,n,l,k,m)=zz p(i,n,l,m,k)=zz p(i,n,m,k,l)=zz p(i,n,m,l,k)=zz 51 continue c 21 continue c return end c c_____________________________________________________________________________ c SUBROUTINE ESUM(par) c c c Sum of electrostatic interaction energy over pairs of sites and hard c sphere repulsion c implicit real*8 (a-h,o-z) common /blk0/q0(50),q1(50,3),q2(50,5),q3(50,7),q4(50,9), . vdwrad(50),x0(50,3),d0(50,3),t0(50,3,3),u0(50,3,3,3) common /blk1/del(3,3),x(50,3),q(50), . d(50,3),t(50,3,3),u(50,3,3,3),radius, . dtot(3),ttot(3,3),utot(3,3,3),xa(3),xab(3), . c(3,3),an(3),pi,conv,na,nb,esint(5,5),reptrm, . ilevel,ifit common /hexad/p(50,3,3,3,3),p0(50,3,3,3,3) common /cdum/ndummy,cdummy(10,4) common /ttblk/tt(6561) common /blkder/iflag equivalence (xx,xab(1)),(yy,xab(2)),(zz,xab(3)) real*8 par(6),r(3) integer*2 ia,ib,ig,id,ie,ip,ic,if,i,j,index,k,l,ii,ij,ik,il C one3=1.d0/3.d0 one9=1.d0/9.d0 one15=1.d0/15.d0 one45=1.d0/45.d0 one105=1.d0/105.d0 one225=1.d0/225.d0 one315=1.d0/315.d0 on1575=1.d0/1575.d0 o11025=1.d0/11025.d0 c reptrm=0.d0 do 241 i=1,5 do 241 j=1,5 esint(i,j)=0.d0 241 continue c c...'dirty' reset of overflowed or underflowed R in the A to B separation c vector and normalisation of Theta and Psi c par(1)=dabs(par(1)) if(ifit.eq.1.and.iflag.le.2.and.par(1).gt.radius) * par(1)=radius if(ifit.eq.1.and.par(1).lt.1.d0)par(1)=5.d0 CALL NORMAL(1,par(2)) CALL NORMAL(2,par(3)) r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) c c...loop over all atoms of molecule A c do 12 i=1,na do 13 ii=1,3 xa(ii)=x(i,ii) 13 continue c c...for every atom of A loop over all atoms of B c do 14 j=na+1, na+nb vcc=0.0d0 vcd=0.0d0 vcq=0.d0 vdd=0.d0 vdq=0.0d0 vco=0.0d0 vqq=0.0d0 vdo=0.0d0 vch=0.0d0 vqo=0.d0 vdh=0.d0 voo=0.d0 vqh=0.d0 voh=0.d0 vhh=0.d0 sepsq=0.0d0 do 15 k=1,3 term=r(k)+x(j,k)-xa(k) sepsq=sepsq+term*term xab(k)=term 15 continue c sep=dsqrt(sepsq) sep4=sepsq*sepsq xx2=xx**2 yy2=yy**2 zz2=zz**2 xx4=xx2*xx2 yy4=yy2*yy2 zz4=zz2*zz2 xxzz=xx*zz xxyy=xx*yy yyzz=yy*zz c c c...T, charge-charge c vcc=vcc+q(i)*q(j)/sep c c c...T.a, charge-dipole c const=1.d0/sep**3 do 16 k=1,3 t1=-xab(k)*const vcd=vcd+t1*(q(i)*d(j,k)-q(j)*d(i,k)) 16 continue c c c...T.ab, charge-quadrupole and dipole-dipole c const=1.d0/sep**5 do 17 k=1,3 do 17 l=1,3 t2=(3.d0*xab(k)*xab(l)-sepsq*del(k,l)) . *const vcq=vcq+t2*(q(i)*t(j,k,l)+q(j)*t(i,k,l)) vdd=vdd-t2*d(i,k)*d(j,l) 17 continue vcq=vcq*one3 c c c...T.abg, dipole-quadrupole + charge-octopole c const=-3.d0/sep**7 do 18 ii=1,3 do 18 ij=1,3 do 18 ik=1,3 terma=xab(ii)*del(ij,ik)+xab(ij)*del(ii,ik) . +xab(ik)*del(ii,ij) t3=(5.0d0*xab(ii)*xab(ij)*xab(ik) . -sepsq*terma)*const vdq=vdq+t3*(-d(i,ii)*t(j,ij,ik) . +d(j,ii)*t(i,ij,ik)) if(ilevel.eq.1)goto 18 vco=vco+t3*(q(i)*u(j,ii,ij,ik) . -q(j)*u(i,ii,ij,ik)) 18 continue vdq=vdq*one3 vco=vco*one15 c c c...T.abgd, quadrupole-quadrupole + dipole-octopole + charge-hexadecapole c c (some increase in speed is obtained by evaluating the nonequivalent c components of T.abgd prior to summing the energy terms) c c const=3.d0/sep**9 const1=5.d0*const ssq3=3.d0*sepsq ssq5=5.d0*sepsq ssq30=30.d0*sepsq ssf3=3.d0*sep4 xs7=7.d0*xx2 ys7=7.d0*yy2 zs7=7.d0*zz2 tt(1)= const*(35.d0*xx4-ssq30*xx2+ssf3) tt(16)=const*(35.d0*yy4-ssq30*yy2+ssf3) tt(81)=const*(35.d0*zz4-ssq30*zz2+ssf3) c tt(2)= const1*xxyy*(xs7-ssq3) tt(24)=const1*yyzz*(ys7-ssq3) tt(8)= const1*xxyy*(ys7-ssq3) tt(3)= const1*xxzz*(xs7-ssq3) tt(27)=const1*xxzz*(zs7-ssq3) tt(54)=const1*yyzz*(zs7-ssq3) c tt(4)= const*(5.d0*xx2*(ys7-sepsq)-ssq5*yy2+sep4) tt(36)=const*(5.d0*yy2*(zs7-sepsq)-ssq5*zz2+sep4) tt(9)= const*(5.d0*xx2*(zs7-sepsq)-ssq5*zz2+sep4) c tt(18)=const1*xxyy*(zs7-sepsq) tt(12)=const1*xxzz*(ys7-sepsq) tt(6)= const1*yyzz*(xs7-sepsq) c do 19 ii=1,3 do 19 ij=1,3 do 19 ik=1,3 do 19 il=1,3 index=ii*ij*ik*il vqq=vqq+tt(index)*t(i,ii,ij)*t(j,ik,il) goto(19,100,101,100,100,100),ilevel 100 vch=vch+tt(index)*(q(i)*p(j,ii,ij,ik,il) . +q(j)*p(i,ii,ij,ik,il)) vdo=vdo+tt(index)*(d(i,ii)*u(j,ij,ik,il) . +d(j,ii)*u(i,ij,ik,il)) goto 19 101 vdo=vdo+tt(index)*(d(i,ii)*u(j,ij,ik,il) . +d(j,ii)*u(i,ij,ik,il)) 19 continue vqq=vqq*one9 vdo=-vdo*one15 vch=vch*one105 c c c...T.abgde, quadrupole-octopole + dipole-hexadecapole c if(ilevel.le.2)goto 71 c xxyyzz=xxyy*zz const=-15.d0/sep**11 const1=3.d0*const const2=21.d0*const ssq7=7.d0*sepsq ssq14=14.d0*sepsq ssq70=70.d0*sepsq ssf15=15.d0*sep4 xs3=3.d0*xx2 ys3=3.d0*yy2 zs3=3.d0*zz2 xs63=63.d0*xx2 ys63=63.d0*yy2 c tt(1)= const*xx*(63.d0*xx4-ssq70*xx2+ssf15) tt(32)= const*yy*(63.d0*yy4-ssq70*yy2+ssf15) tt(243)=const*zz*(63.d0*zz4-ssq70*zz2+ssf15) c term=const1*(21.d0*xx4-ssq14*xx2+sep4) tt(2)= yy*term tt(3)= zz*term term=const1*(21.d0*yy4-ssq14*yy2+sep4) tt(16)=xx*term tt(48)=zz*term term=const1*(21.d0*zz4-ssq14*zz2+sep4) tt(81)=xx*term tt(162)=yy*term c tt(4)= const*xx*(ssf3-ssq7*(xx2+ys3)+xs63*yy2) tt(9)= const*xx*(ssf3-ssq7*(xx2+zs3)+xs63*zz2) tt(8)= const*yy*(ssf3-ssq7*(xs3+yy2)+xs63*yy2) tt(72)=const*yy*(ssf3-ssq7*(yy2+zs3)+ys63*zz2) tt(27)=const*zz*(ssf3-ssq7*(xs3+zz2)+xs63*zz2) tt(108)=const*zz*(ssf3-ssq7*(ys3+zz2)+ys63*zz2) c tt(6)= const2*xxyyzz*(xs3-sepsq) tt(24)=const2*xxyyzz*(ys3-sepsq) tt(54)=const2*xxyyzz*(zs3-sepsq) c tt(12)=const*zz*(sep4-ssq7*(xx2+yy2)+xs63*yy2) tt(18)=const*yy*(sep4-ssq7*(xx2+zz2)+xs63*zz2) tt(36)=const*xx*(sep4-ssq7*(yy2+zz2)+ys63*zz2) c do 50 ia=1,3 do 50 ib=1,3 do 50 ig=1,3 do 50 id=1,3 do 50 ie=1,3 index=ia*ib*ig*id*ie vqo=vqo+tt(index)*(t(i,ia,ib)*u(j,ig,id,ie) . -t(j,ia,ib)*u(i,ig,id,ie)) if(ilevel.lt.4)goto 50 vdh=vdh+tt(index)*(-d(i,ia)*p(j,ib,ig,id,ie) . +d(j,ia)*p(i,ib,ig,id,ie)) 50 continue vqo=vqo*one45 vdh=vdh*one105 c c c...T.abgdep, octopole-octopole + quadrupole-hexadecapole c const=1.d0/sep**13 const1=315.d0*const const2=945.d0*const ssq18=18.d0*sepsq ssq30=30.d0*sepsq ssq945=945.d0*sepsq ssql=14175.d0*sepsq ssf5=5.d0*sep4 ssf315=315.d0*sep4 ssfl=4725.d0*sep4 sep6=sep4*sepsq sssm45=-45.d0*sep6 sss225=225.d0*sep6 xs33=33.d0*xx2 ys33=33.d0*yy2 xf33=33.d0*xx4 yf33=33.d0*yy4 zf33=33.d0*zz4 xsl=10395.d0*xx2 ysl=10395.d0*yy2 zsl=10395.d0*zz2 xs2=2.d0*xx2 ys2=2.d0*yy2 zs2=2.d0*zz2 xs6=6.d0*xx2 ys6=6.d0*yy2 zs6=6.d0*zz2 c tt(1)= const*(xsl*xx4-ssql*xx4+ssfl*xx2-sss225) tt(64)= const*(ysl*yy4-ssql*yy4+ssfl*yy2-sss225) tt(729)= const*(zsl*zz4-ssql*zz4+ssfl*zz2-sss225) c term= const1*(xf33-ssq30*xx2+ssf5) tt(2)=term*xxyy tt(3)=term*xxzz term= const1*(yf33-ssq30*yy2+ssf5) tt(32)=term*xxyy tt(96)=term*yyzz term= const1*(zf33-ssq30*zz2+ssf5) tt(243)=term*xxzz tt(486)=term*yyzz c tt(4)= const*(sssm45+ssf315*(xs2+yy2) . -ssq945*xx2*(xx2+ys6)+ysl*xx4) tt(9)= const*(sssm45+ssf315*(xs2+zz2) . -ssq945*xx2*(xx2+zs6)+zsl*xx4) tt(16)=const*(sssm45+ssf315*(xx2+ys2) . -ssq945*yy2*(xs6+yy2)+xsl*yy4) tt(144)=const*(sssm45+ssf315*(ys2+zz2) . -ssq945*yy2*(yy2+zs6)+zsl*yy4) tt(81)=const*(sssm45+ssf315*(xx2+zs2) . -ssq945*zz2*(xs6+zz2)+xsl*zz4) tt(324)=const*(sssm45+ssf315*(yy2+zs2) . -ssq945*zz2*(ys6+zz2)+ysl*zz4) c tt(8)= const2*xxyy*(sep4-ssq3*(xx2+yy2)+11.d0*xx2*yy2) tt(27)=const2*xxzz*(sep4-ssq3*(xx2+zz2)+11.d0*xx2*zz2) tt(216)=const2*yyzz*(sep4-ssq3*(yy2+zz2)+11.d0*yy2*zz2) c tt(6)= const1*yyzz*(xf33-ssq18*xx2+sep4) tt(48)=const1*xxzz*(yf33-ssq18*yy2+sep4) tt(162)=const1*xxyy*(zf33-ssq18*zz2+sep4) c tt(12)= const1*xxzz*(sep4-ssq3*(xx2+ys3)+xs33*yy2) tt(18)= const1*xxyy*(sep4-ssq3*(xx2+zs3)+xs33*zz2) tt(24)= const1*yyzz*(sep4-ssq3*(xs3+yy2)+xs33*yy2) tt(72)= const1*xxyy*(sep4-ssq3*(yy2+zs3)+ys33*zz2) tt(54)= const1*yyzz*(sep4-ssq3*(xs3+zz2)+xs33*zz2) tt(108)= const1*xxzz*(sep4-ssq3*(ys3+zz2)+ys33*zz2) c tt(36)=const*(90.d0*sep6-ssq945*(xx2*(yy2+zz2) . +yy2*zz2)+xsl*yy2*zz2) c do 51 ia=1,3 do 51 ib=1,3 do 51 ig=1,3 do 51 id=1,3 do 51 ie=1,3 do 51 ip=1,3 index=ia*ib*ig*id*ie*ip voo=voo+tt(index) . *u(i,ia,ib,ig)*u(j,id,ie,ip) if(ilevel.lt.4)goto 51 vqh=vqh+tt(index) . *(t(i,ia,ib)*p(j,ig,id,ie,ip) . +t(j,ia,ib)*p(i,ig,id,ie,ip)) 51 continue voo=-voo*one225 vqh=vqh*one315 c c c...T.abgdepc, octopole-hexadecapole c if(ilevel.lt.5)goto 71 xyz429=429.d0*xxyyzz*xxyyzz xx6c=429.d0*xx2*xx4 yy6c=429.d0*yy2*yy4 zz6c=429.d0*zz2*zz4 ssq33=33.d0*sepsq ssq110=110.d0*sepsq ssq495=495.d0*sepsq ssq693=693.d0*sepsq ssf9=9.d0*sep4 ssf39=39.d0*sep4 ssf135=135.d0*sep4 sss5=5.d0*sep6 sss8=8.d0*sep6 sss35=35.d0*sep6 const=315.d0/sep**15 constm=-const const1=3.d0*const const2=-const1 xs10=10.d0*xx2 ys10=10.d0*yy2 zs10=10.d0*zz2 xs143=143.d0*xx2 ys143=143.d0*yy2 zs143=143.d0*zz2 xs429=429.d0*xx2 ys429=429.d0*yy2 zs429=429.d0*zz2 c t_aaaaaaa: tt( 1)=const*xx*(sss35-ssf315*xx2+ssq693*xx4-xx6c) xxxxxxx tt( 128)=const*yy*(sss35-ssf315*yy2+ssq693*yy4-yy6c) yyyyyyy tt(2187)=const*zz*(sss35-ssf315*zz2+ssq693*zz4-zz6c) zzzzzzz c t_aaaaaab: tt( 2)=const*yy*(sss5-ssf135*xx2+ssq495*xx4-xx6c) xxxxxxy tt( 3)=const*zz*(sss5-ssf135*xx2+ssq495*xx4-xx6c) xxxxxxz tt( 64)=const*xx*(sss5-ssf135*yy2+ssq495*yy4-yy6c) yyyyyyx tt( 192)=const*zz*(sss5-ssf135*yy2+ssq495*yy4-yy6c) yyyyyyz tt( 729)=const*xx*(sss5-ssf135*zz2+ssq495*zz4-zz6c) zzzzzzx tt(1458)=const*yy*(sss5-ssf135*zz2+ssq495*zz4-zz6c) zzzzzzy c t_aaaaabb: tt( 4)=const*xx*(sss5-ssf15*(xs2+ys3) xxxxxyy . +ssq33*xx2*(xx2+ys10)-ys429*xx4) tt( 9)=const*xx*(sss5-ssf15*(xs2+zs3) xxxxxzz . +ssq33*xx2*(xx2+zs10)-zs429*xx4) tt( 32)=const*yy*(sss5-ssf15*(xs3+ys2) yyyyyxx . +ssq33*yy2*(xs10+yy2)-xs429*yy4) tt(288)=const*yy*(sss5-ssf15*(ys2+zs3) yyyyyzz . +ssq33*yy2*(yy2+zs10)-zs429*yy4) tt(243)=const*zz*(sss5-ssf15*(xs3+zs2) zzzzzxx . +ssq33*zz2*(xs10+zz2)-xs429*zz4) tt(972)=const*zz*(sss5-ssf15*(ys3+zs2) zzzzzyy . +ssq33*zz2*(ys10+zz2)-ys429*zz4) c t_aaaabbb: tt( 8)=const1*yy*(sep6-ssf3*(xs6+yy2) xxxxyyy . +ssq33*xx2*(xx2+ys2)-ys143*xx4) tt( 27)=const1*zz*(sep6-ssf3*(xs6+zz2) xxxxzzz . +ssq33*xx2*(xx2+zs2)-zs143*xx4) tt( 16)=const1*xx*(sep6-ssf3*(xx2+ys6) yyyyxxx . +ssq33*yy2*(xs2+yy2)-xs143*yy4) tt(432)=const1*zz*(sep6-ssf3*(ys6+zz2) yyyyzzz . +ssq33*yy2*(yy2+zs2)-zs143*yy4) tt( 81)=const1*xx*(sep6-ssf3*(xx2+zs6) zzzzxxx . +ssq33*zz2*(xs2+zz2)-xs143*zz4) tt(648)=const1*yy*(sep6-ssf3*(yy2+zs6) zzzzyyy . +ssq33*zz2*(ys2+zz2)-ys143*zz4) c t_aaaaabg: tt( 6)=const2*xxyyzz*(ssf15-ssq110*xx2+143*xx4) xxxxxyz tt( 96)=const2*xxyyzz*(ssf15-ssq110*yy2+143*yy4) yyyyyxz tt(486)=const2*xxyyzz*(ssf15-ssq110*zz2+143*zz4) zzzzzxy c t_aaaabbg: tt( 12)=const*zz*(sep6-ssf9*(xs2+yy2) xxxxyyz . +ssq33*xx2*(xx2+ys6)-ys429*xx4) tt( 18)=const*yy*(sep6-ssf9*(xs2+zz2) xxxxzzy . +ssq33*xx2*(xx2+zs6)-zs429*xx4) tt( 48)=const*zz*(sep6-ssf9*(xx2+ys2) yyyyxxz . +ssq33*yy2*(xs6+yy2)-xs429*yy4) tt(144)=const*xx*(sep6-ssf9*(ys2+zz2) yyyyzzx . +ssq33*yy2*(yy2+zs6)-zs429*yy4) tt(162)=const*yy*(sep6-ssf9*(xx2+zs2) zzzzxxy . +ssq33*zz2*(xs6+zz2)-xs429*zz4) tt(324)=const*xx*(sep6-ssf9*(yy2+zs2) zzzzyyx . +ssq33*zz2*(ys6+zz2)-ys429*zz4) c t_aaabbbg: tt( 24)=const2*xxyyzz*(ssf9-ssq33*(xx2+yy2)+xs143*yy2) xxxyyyz tt( 54)=const2*xxyyzz*(ssf9-ssq33*(xx2+zz2)+xs143*zz2) xxxzzzy tt(216)=const2*xxyyzz*(ssf9-ssq33*(yy2+zz2)+ys143*zz2) yyyzzzx c t_aaabbgg: tt( 36)=constm*xx*(sss8-ssf39*xx2+ssq33*(xx4-ys3*zz2)+xyz429) xxxyyzz tt( 72)=constm*yy*(sss8-ssf39*yy2-ssq33*(xs3*zz2-yy4)+xyz429) yyyxxzz tt(108)=constm*zz*(sss8-ssf39*zz2-ssq33*(xs3*yy2-zz4)+xyz429) zzzxxyy c do 80 ia=1,3 do 80 ib=1,3 do 80 ig=1,3 do 80 id=1,3 do 80 ie=1,3 do 80 ip=1,3 do 80 ic=1,3 index=ia*ib*ig*id*ie*ip*ic voh=voh+tt(index) . *(u(i,ia,ib,ig)*p(j,id,ie,ip,ic) . -u(j,ia,ib,ig)*p(i,id,ie,ip,ic)) 80 continue voh=-voh*on1575 c c c...T.abgdepcf, hexadecapole-hexadecapole c if(ilevel.ne.6)goto 71 const=315.d0/sep**17 constm=-const const1=3.d0*const con1m=-const1 const2=-9.d0*const const3=-45.d0*const sep8=sep4*sep4 sse35=35.d0*sep8 sse5=5.d0*sep8 sse8=8.d0*sep8 sss45=45.d0*sep6 sss15=15.d0*sep6 sss18=18.d0*sep6 sss90=90.d0*sep6 sss3=3.d0*sep6 sss9=9.d0*sep6 sssl=1260.d0*sep6 ssf11=11.d0*sep4 ssf33=33.d0*sep4 ssf55=55.d0*sep4 ssf99=99.d0*sep4 ssf165=165.d0*sep4 ssf385=385.d0*sep4 ssf495=495.d0*sep4 ssfl=6930.d0*sep4 ssq143=143.d0*sepsq ssq429=429.d0*sepsq ssq858=858.d0*sepsq ssqt1=1001.d0*sepsq ssql=12012.d0*sepsq xx6=xx2*xx4 yy6=yy2*yy4 zz6=zz2*zz4 xf3=3.d0*xx4 yf3=3.d0*yy4 zf3=3.d0*zz4 xf6=6.d0*xx4 yf6=6.d0*yy4 zf6=6.d0*zz4 xf2145=2145.d0*xx4 yf2145=2145.d0*yy4 zf2145=2145.d0*zz4 xf6435=6435.d0*xx4 yf6435=6435.d0*yy4 zf6435=6435.d0*zz4 xxyys=xxyy*xxyy xxzzs=xxzz*xxzz yyzzs=yyzz*yyzz xs15=15.d0*xx2 ys15=15.d0*yy2 zs15=15.d0*zz2 xs715=715.d0*xx2 ys715=715.d0*yy2 zs715=715.d0*zz2 xss715=715.d0*xx6 yss715=715.d0*yy6 zss715=715.d0*zz6 c tt( 1)=const*(sse35-sssl*xx2+ssfl*xx4 . -ssql*xx6+xf6435*xx4) tt( 256)=const*(sse35-sssl*yy2+ssfl*yy4 . -ssql*yy6+yf6435*yy4) tt(6561)=const*(sse35-sssl*zz2+ssfl*zz4 . -ssql*zz6+zf6435*zz4) c tt( 2)=const2*xxyy*(sss35-ssf385*xx2+ssqt1*xx4-xss715) tt( 3)=const2*xxzz*(sss35-ssf385*xx2+ssqt1*xx4-xss715) tt( 128)=const2*xxyy*(sss35-ssf385*yy2+ssqt1*yy4-yss715) tt( 384)=const2*yyzz*(sss35-ssf385*yy2+ssqt1*yy4-yss715) tt(2187)=const2*xxzz*(sss35-ssf385*zz2+ssqt1*zz4-zss715) tt(4374)=const2*yyzz*(sss35-ssf385*zz2+ssqt1*zz4-zss715) c tt( 4)=const*(sse5-sss45*(xs3+yy2)+ssf495*xx2*(xx2+ys3) . -ssq429*xx4*(xx2+ys15)+xf6435*xxyys) tt( 9)=const*(sse5-sss45*(xs3+zz2)+ssf495*xx2*(xx2+zs3) . -ssq429*xx4*(xx2+zs15)+xf6435*xxzzs) tt( 64)=const*(sse5-sss45*(xx2+ys3)+ssf495*yy2*(xs3+yy2) . -ssq429*yy4*(xs15+yy2)+yf6435*xxyys) tt( 576)=const*(sse5-sss45*(ys3+zz2)+ssf495*yy2*(yy2+zs3) . -ssq429*yy4*(yy2+zs15)+yf6435*yyzzs) tt( 729)=const*(sse5-sss45*(xx2+zs3)+ssf495*zz2*(xs3+zz2) . -ssq429*zz4*(xs15+zz2)+zf6435*xxzzs) tt(2916)=const*(sse5-sss45*(yy2+zs3)+ssf495*zz2*(ys3+zz2) . -ssq429*zz4*(ys15+zz2)+zf6435*yyzzs) c tt( 8)=con1m*xxyy*(sss45-ssf165*(xs2+yy2) . +ssq143*xx2*(xs3+ys10)-xf2145*yy2) tt( 27)=con1m*xxzz*(sss45-ssf165*(xs2+zz2) . +ssq143*xx2*(xs3+zs10)-xf2145*zz2) tt( 32)=con1m*xxyy*(sss45-ssf165*(xx2+ys2) . +ssq143*yy2*(xs10+ys3)-yf2145*xx2) tt( 864)=con1m*yyzz*(sss45-ssf165*(ys2+zz2) . +ssq143*yy2*(ys3+zs10)-yf2145*zz2) tt( 243)=con1m*xxzz*(sss45-ssf165*(xx2+zs2) . +ssq143*zz2*(xs10+zs3)-zf2145*xx2) tt(1944)=con1m*yyzz*(sss45-ssf165*(yy2+zs2) . +ssq143*zz2*(ys10+zs3)-zf2145*yy2) c tt( 16)=const1*(sep8-sss18*(xx2+yy2)+ssf33*(xx4 . +12.d0*xxyys+yy4)-ssq858*xxyys*(xx2+yy2)+xf2145*yy4) tt( 81)=const1*(sep8-sss18*(xx2+zz2)+ssf33*(xx4 . +12.d0*xxzzs+zz4)-ssq858*xxzzs*(xx2+zz2)+xf2145*zz4) tt(1296)=const1*(sep8-sss18*(yy2+zz2)+ssf33*(yy4 . +12.d0*yyzzs+zz4)-ssq858*yyzzs*(yy2+zz2)+yf2145*zz4) tt( 6)=const3*yyzz*(sep6-ssf33*xx2+ssq143*xx4-143.d0*xx6) tt( 192)=const3*xxzz*(sep6-ssf33*yy2+ssq143*yy4-143.d0*yy6) tt(1458)=const3*xxyy*(sep6-ssf33*zz2+ssq143*zz4-143.d0*zz6) c tt( 12)=con1m*xxzz*(sss15-ssf55*(xs2+ys3) . +ssq143*xx2*(xx2+ys10)-xf2145*yy2) tt( 18)=con1m*xxyy*(sss15-ssf55*(xs2+zs3) . +ssq143*xx2*(xx2+zs10)-xf2145*zz2) tt( 96)=con1m*yyzz*(sss15-ssf55*(xs3+ys2) . +ssq143*yy2*(xs10+yy2)-yf2145*xx2) tt( 288)=con1m*xxyy*(sss15-ssf55*(ys2+zs3) . +ssq143*yy2*(yy2+zs10)-yf2145*zz2) tt( 486)=con1m*yyzz*(sss15-ssf55*(xs3+zs2) . +ssq143*zz2*(xs10+zz2)-zf2145*xx2) tt( 972)=con1m*xxzz*(sss15-ssf55*(ys3+zs2) . +ssq143*zz2*(ys10+zz2)-zf2145*yy2) c tt( 24)=const2*yyzz*(sss3-ssf11*(xs6+yy2) . +ssq143*xx2*(xx2+ys2)-ys715*xx4) tt( 54)=const2*yyzz*(sss3-ssf11*(xs6+zz2) . +ssq143*xx2*(xx2+zs2)-zs715*xx4) tt( 48)=const2*xxzz*(sss3-ssf11*(xx2+ys6) . +ssq143*yy2*(xs2+yy2)-xs715*yy4) tt( 432)=const2*xxzz*(sss3-ssf11*(ys6+zz2) . +ssq143*yy2*(yy2+zs2)-zs715*yy4) tt( 648)=const2*xxyy*(sss3-ssf11*(yy2+zs6) . +ssq143*zz2*(ys2+zz2)-ys715*zz4) tt( 162)=const2*xxyy*(sss3-ssf11*(xx2+zs6) . +ssq143*zz2*(xs2+zz2)-xs715*zz4) c tt( 36)=constm*(sse8-sss9*(21.d0*xx2+11.d0*yy2) . +ssf99*(xf6+27.d0*xxyys+yy4)-ssq429*xx2 . *(xx4+6.d0*xxyys+yf6)-xf6435*yyzzs) tt( 144)=constm*(sse8-sss9*(11.d0*xx2+21.d0*yy2) . +ssf99*(xx4+27.d0*xxyys+yf6)-ssq429*yy2 . *(xf6+6.d0*xxyys+yy4)-yf6435*xxzzs) tt( 324)=constm*(sse8-sss9*(11.d0*xx2+21.d0*zz2) . +ssf99*(xx4+27.d0*xxzzs+zf6)-ssq429*zz2 . *(xf6+6.d0*xxzzs+zz4)-zf6435*xxyys) c tt( 72)=const1*xxyy*(sss90-ssf495*(xx2+yy2)+ssq143*(xf3 . +20.d0*xxyys+yf3)-2145.d0*xxyys*(xx2+yy2)) tt( 108)=const1*xxzz*(sss90-ssf495*(xx2+zz2)+ssq143*(xf3 . +20.d0*xxzzs+zf3)-2145.d0*xxzzs*(xx2+zz2)) tt( 216)=const1*yyzz*(sss90-ssf495*(yy2+zz2)+ssq143*(yf3 . +20.d0*yyzzs+zf3)-2145.d0*yyzzs*(yy2+zz2)) c do 85 ia=1,3 do 85 ib=1,3 do 85 ig=1,3 do 85 id=1,3 do 85 ie=1,3 do 85 ip=1,3 do 85 ic=1,3 do 85 if=1,3 index=ia*ib*ig*id*ie*ip*ic*if vhh=vhh+tt(index) . *p(i,ia,ib,ig,id)*p(j,ie,ip,ic,if) 85 continue vhh=vhh*o11025 c 71 esint(1,1)=esint(1,1)+vcc esint(2,1)=esint(2,1)+vcd esint(2,2)=esint(2,2)+vdd esint(3,1)=esint(3,1)+vcq esint(3,2)=esint(3,2)+vdq esint(3,3)=esint(3,3)+vqq esint(4,1)=esint(4,1)+vco esint(4,2)=esint(4,2)+vdo esint(4,3)=esint(4,3)+vqo esint(4,4)=esint(4,4)+voo esint(5,1)=esint(5,1)+vch esint(5,2)=esint(5,2)+vdh esint(5,3)=esint(5,3)+vqh esint(5,4)=esint(5,4)+voh esint(5,5)=esint(5,5)+vhh c c c...repulsive term c sigij=vdwrad(i)+vdwrad(j) if(sigij.eq.0.d0)goto 14 if(sep.gt.2.d0*sigij)goto 14 c if(sep.ge.(sigij-0.6d0))then arg=(sigij-sep-0.6d0)/0.1d0 reptrm=reptrm+25.d0*(1.d0+dtanh(arg)) c else if(sep.lt.(sigij-0.6d0))then c reptrm=reptrm+250.d0*(sigij-0.6d0-sep+1.d0)**10-225.d0 c endif c 14 continue 12 continue c c...repulsive contribution from dummy centres c if(ndummy.eq.0)goto 49 do 61 i=1,ndummy do 60 j=na+1,na+nb sepsq=0.d0 do 62 k=1,3 sepsq=sepsq+(x(j,k)+r(k)-cdummy(i,k))**2 62 continue sep=dsqrt(sepsq) sigij=cdummy(i,4)+vdwrad(j) if(sep.gt.2.d0*sigij)goto 60 if(sep.ge.(sigij-0.6d0))then arg=(sigij-sep-0.6d0)/0.1d0 reptrm=reptrm+25.d0*(1.d0+dtanh(arg)) else if(sep.lt.(sigij-0.6d0))then reptrm=reptrm+25.d0+250.d0*(sigij-0.5d0)/(sep+0.1d0) endif 60 continue 61 continue c 49 return end c_____________________________________________________________________________ c SUBROUTINE DFPMIN(P,N,FTOL,ITER,FRET) C C Given a starting point P that is a vector of length N, the Broyden- C Fletcher-Goldfarb-Shanno varian of Davidon-Fletcher-Powell minimisation C is performed on a functin FUNC, using its gradient as calculated by a C routine DFUNC. The convergence requirement on the function value is C input as FTOL. Returned quantities are P (the location of the minimum), C ITER (the number of iterations that were performed), and FRET (the minimum C value of the function). The routine LINMIN is called to perform line C minimisations. C C Minimisation routines: DFPMIN, LINMIN, MNBRAK, BRENT and F1DIM C come from Numerical Recipes. Modified or added lines are indicated by a C 'ZK' in the comment field. C C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6,ITMAX=100,EPS=1.d-10) common /mon/imonit,filmon ZK DIMENSION P(6),HESSIN(NMAX,NMAX),XI(NMAX),G(NMAX),DG(NMAX), ZK * HDG(NMAX) character*25 filmon ZK pi=3.1415926535897932d0 ZK conv=pi/180d0 ZK bohr=0.5291771d0 ZK C C...calculate starting function value and gradient C FP=ESINTE(P) ZK CALL ESINTD(P,G,FP) ZK write(*,30)0,fp,fp,p(1)*bohr,(p(i)/conv,i=2,n), ZK * g(1)*bohr,(g(i)/conv,i=2,n) ZK if(imonit.eq.3)then ZK open(4,file=filmon,access='APPEND') ZK write(4,30)0,fp,fp,p(1)*bohr,(p(i)/conv,i=2,n), ZK * g(1)*bohr,(g(i)/conv,i=2,n) ZK close(4) ZK endif ZK C C...initialise the inverse Hessian to the unit matrix and set initial C line direction XI C DO 12 I=1,N DO 11 J=1,N HESSIN(I,J)=0.d0 11 CONTINUE HESSIN(I,I)=1.d0 XI(I)=-G(I) 12 CONTINUE C C...main loop over the iterations C DO 24 ITER=1,ITMAX CALL LINMIN(P,XI,N,FRET) if(imonit.eq.3)then open(4,file=filmon,access='APPEND') ZK write(4,30)iter,fcalc,fret,p(1)*bohr, ZK * (p(i)/conv,i=2,n),g(1)*bohr,(g(i)/conv,i=2,n) ZK close(4) ZK endif ZK write(*,30)iter,fcalc,fret,p(1)*bohr,(p(i)/conv,i=2,n), ZK * g(1)*bohr,(g(i)/conv,i=2,n) ZK 30 format(' iteration no.',i2,' E.calc = ',f16.10, ZK * ' E.pred =',f16.10 /1x,6f12.6/1x,6f12.6) ZK C C...return if converged to a minimum IF(2.d0*DABS(FRET-FP).LE.FTOL*(DABS(FRET)+DABS(FP)+EPS)) * RETURN C C...save old function value and the old gradient FP=FRET DO 13 I=1,N DG(I)=G(I) 13 CONTINUE C C...get new function value and gradient FRET=ESINTE(P) ZK CALL ESINTD(P,G,FRET) ZK fcalc=fret ZK C C...compute difference of gradients DO 14 I=1,N DG(I)=G(I)-DG(I) 14 CONTINUE C C...and difference times the current matrix DO 16 I=1,N HDG(I)=0.d0 DO 15 J=1,N HDG(I)=HDG(I)+HESSIN(I,J)*DG(J) 15 CONTINUE 16 CONTINUE C C...calculate dot products for the denominators FAC=0.d0 FAE=0.d0 DO 17 I=1,N FAC=FAC+DG(I)*XI(I) FAE=FAE+DG(I)*HDG(I) 17 CONTINUE C C...and make the denominators multiplicative FAC=1.d0/FAC FAD=1.d0/FAE C C...the vector which makes BFGS different from DFP DO 18 I=1,N DG(I)=FAC*XI(I)-FAD*HDG(I) 18 CONTINUE C C...the BFGS updating formula DO 21 I=1,N DO 19 J=1,N HESSIN(I,J)=HESSIN(I,J)+FAC*XI(I)*XI(J) * -FAD*HDG(I)*HDG(J)+FAE*DG(I)*DG(J) 19 CONTINUE 21 CONTINUE C C...now calculate the next direction to go DO 23 I=1,N XI(I)=0.d0 DO 22 J=1,N XI(I)=XI(I)-HESSIN(I,J)*G(J) 22 CONTINUE 23 CONTINUE C C...and go back to another iteration 24 CONTINUE C WRITE(*,'(1x/a/)')'*********** Too many iterations in DFPMIN' open(4,file=filmon,access='APPEND') ZK WRITE(4,'(1x/a)')'*********** Too many iterations in DFPMIN' ZK close(4) ZK C RETURN END C C_____________________________________________________________________________ C SUBROUTINE NORMAL(I,ANGLE) C C Normalisation of an angle (in radians) to range -PI to PI for I=1 and C 0 to 2*PI for I=2 C IMPLICIT REAL*8 (A-H,O-Z) DATA PI,TPI/3.141592654D0,6.283185307D0/ C IF(I.EQ.1)THEN A=DABS(ANGLE) IF(A.LE.PI)RETURN A=A-DINT(A/TPI)*TPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.-PI)ANGLE=ANGLE+TPI IF(ANGLE.GT.PI)ANGLE=ANGLE-TPI RETURN ELSE IF(I.EQ.2)THEN A=DABS(ANGLE) A=A-DINT(A/TPI)*TPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.0.D0)ANGLE=ANGLE+TPI RETURN ENDIF C RETURN END C C_____________________________________________________________________________ C SUBROUTINE LINMIN(P,XI,N,FRET) C C Given an N dimensional point P and an N dimensional direction XI, LINMIN C moves and resets P to where the function FUNC(P) takes on a minimum C along the direction XI from P, and replaces XI by the actual vector C displacement that P was moved. Also returns FRET the value of FUNC at C the returned location P. This is actually all accomplished by calling C the routines MNBRAK and BRENT. C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6,TOL=3.d-4) EXTERNAL F1DIM DIMENSION P(6),XI(6) ZK COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX) C C...set up the common block C NCOM=N DO 11 J=1,N PCOM(J)=P(J) XICOM(J)=XI(J) 11 CONTINUE C C...initial guess for brackets C AX=0.d0 XX=1.d0 BX=2.d0 CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM) FRET=BRENT(AX,XX,BX,F1DIM,TOL,XMIN) C C...construct the vector results to return C DO 12 J=1,N XI(J)=XMIN*XI(J) P(J)=P(J)+XI(J) 12 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC) C C Given a function FUNC, and given distinct initial points AX and BX, this C routine searches in the downhill direction (defined by the function as C evaluated at the initial points) and returns new points AX, BX, CX which C bracket a minimum of the function. Also returned are the function values C at the three points , FA, FB, and FC. C IMPLICIT REAL*8 (A-H,O-Z) ZK C C...The first parameter is the default ratio by which successive intervals C are magnified; the second is the maximum magnification allowed for a C parabolic-fit step C PARAMETER (GOLD=1.618034d0, GLIMIT=100.d0, TINY=1.d-20) C FA=FUNC(AX) FB=FUNC(BX) C C...Switch roles of A and B so we can go downhill in the direction from A to B C IF(FB.GT.FA)THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF C C...First guess for C C CX=BX+GOLD*(BX-AX) FC=FUNC(CX) C C...'DO WHILE': keep returning here until we bracket C 1 IF(FB.GE.FC)THEN C C...compute U by parabolic extrapolation from A,B,C. TINY is used to prevent C any possible division by zero R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R) * /(2.d0*DSIGN(DMAX1(DABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) C C...got a minimum between B and C, which will exit IF((BX-U)*(U-CX).GT.0.d0)THEN FU=FUNC(U) IF(FU.LT.FC)THEN AX=BX FA=FB BX=U FB=FU GO TO 1 C C...got a minimum between A and U, which will exit ELSE IF(FU.GT.FB)THEN CX=U FC=FU GO TO 1 ENDIF C C...parabolic fit was no use. Use default magnification U=CX+GOLD*(CX-BX) FU=FUNC(U) C C...parabolic fit is between C and its allowed limit ELSE IF((CX-U)*(U-ULIM).GT.0.d0)THEN FU=FUNC(U) IF(FU.LT.FC)THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=FUNC(U) ENDIF C C...limit parabolic U to maximum allowed value ELSE IF((U-ULIM)*(ULIM-CX).GE.0.d0)THEN U=ULIM FU=FUNC(U) C C...reject parabolic U, use default magnification ELSE U=CX+GOLD*(CX-BX) FU=FUNC(U) ENDIF C C...eliminate oldest point and continue AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GO TO 1 ENDIF C RETURN END C C_____________________________________________________________________________ C FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN) C C Given a function F, and given a bracketing triplet of abscissas AX, BX, C CX (such that BX is between AX and CX, and F(BX) is less than both F(AX) C and F(CX)), this routine isolates the minimum to a fractional precision C of about TOL using Brent's method. The abscissa of the minimum is C returned as XMIN, and the minimum function value is returned as BRENT, the C returned function value. C IMPLICIT REAL*8 (A-H,O-Z) ZK C C...maximum allowed number of iterations; golden ratio; and a small number C which protects against trying to achieve fractional accuracy for a minimum C that happens to be exactly zero. C PARAMETER (ITMAX=100,CGOLD=.3819660d0,ZEPS=1.0d-10) C C...A and B must be in ascending order, though the input abscissas need not be C A=DMIN1(AX,CX) B=DMAX1(AX,CX) C C...initialisations C V=BX W=V X=V E=0.d0 FX=F(X) FV=FX FW=FX C C...Main program loop C DO 11 ITER=1,ITMAX XM=0.5d0*(A+B) TOL1=TOL*DABS(X)+ZEPS TOL2=2.d0*TOL1 C C...test for done here IF(DABS(X-XM).LE.(TOL2-.5d0*(B-A))) GOTO 3 C C...construct a trial parabolic fit IF(DABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.d0*(Q-R) IF(Q.GT.0.d0) P=-P Q=DABS(Q) ETEMP=E E=D IF(DABS(P).GE.DABS(.5d0*Q*ETEMP).OR.P.LE.Q*(A-X).OR. * P.GE.Q*(B-X)) GOTO 1 C C...the above conditions determine the acceptability of the parabolic fit. C Here it is OK: take a parabolic step D=P/Q U=X+D IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=DSIGN(TOL1,XM-X) C C...skip over the golden section step GOTO 2 ENDIF C C...we arrive here for a golden section step, which we take into the larger C of the two segments 1 IF(X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF C C...take the golden section step D=CGOLD*E C C...arrive here with D computed either from parabolic fit, or else from golden C section 2 IF(DABS(D).GE.TOL1) THEN U=X+D ELSE U=X+DSIGN(TOL1,D) ENDIF C C...this is the one function evaluation per iteration, and we have to decide C what to do with our function evaluation FU=F(U) C C...housekeeping follows: IF(FU.LE.FX) THEN IF(U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF(U.LT.X) THEN A=U ELSE B=U ENDIF IF(FU.LE.FW .OR. W.EQ.X) THEN V=W FV=FW W=U FW=FU ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V=U FV=FU ENDIF ENDIF C C...done with housekeeping, back for another iteration 11 CONTINUE PAUSE 'Brent exceed maximum iterations.' C C...Arrive here ready to exit with best values C 3 XMIN=X BRENT=FX C RETURN END C C_____________________________________________________________________________ C FUNCTION F1DIM(X) C C This function must accompany LINMIN C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6) COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX) DIMENSION XT(NMAX) C DO 11 J=1,NCOM XT(J)=PCOM(J)+X*XICOM(J) 11 CONTINUE F1DIM=ESINTE(XT) ZK C RETURN END C_____________________________________________________________________________ C_____________________________________________________________________________