c$debug c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c DMAROT - ROTATION OF A DMA ARRAY C C This program reads the data file for the MIN16 program and C allows rotation and translation of the DMA array of molecule B c and outputs the resulting cartesian coordinates and spherical c multipole components. c c Rotation angles are in the same polar coordinate system as c used in MIN16. c c Several rotations can be applied in succession, but it is c probably best to do one translation only, at the end of c rotations. c c Note that shift of DMA origin is up to octupole only. C C c C ver. 28.III.2012 ----- 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 14.12.97: cleaning up (this is the first modification since 2/93) c 31.03.00: translations c 14.04.00: mods to input/output + commenting c 28.03.12: rescued from near oblivion: recompilation and final Cartesians c............................................................................. c c q00 c q10 q11c q11s c q20 q21c q21s q22c q22s c q30 q31c q31s q32c q32s q33c q33s c q40 q41c q41s q42c q42s q43c q43s q44c q44s c............................................................................. c 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) character filnam*25,coment*76,filmon*25,line*80 character*40 label(50) 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,'|'/ * ' | DMAROT - Rotation of multipole array of molecule B', * ' from the .HEX file',T79,'|'/ * ' | defining the DMAs for the A-B system',T79,'|'/ * ' |',76(1H_),'|'/' version 28.III.2012',T64,'Zbigniew KISIEL'/) c 200 write(*,212) * ' Name of the MIN16 data file (with hexadecapoles): ' 212 format(1x/a,$) read(*,'(a)',err=200)filnam open(2,file=filnam,err=200) c 287 write(*,212)' Name of results file: ' read(*,'(a)',err=287)filnam open(3,file=filnam,status='unknown') c c...number of atoms, comment and multipoles for molecule A c read(2,221)na,coment 221 format(i4,a) write(*,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 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) 288 format(' ATOM NUMBER',i3,':') 289 format(a,5f12.6) 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,'(a/)')' UNROTATED DMA OF MOLECULE B:' write(3,221)nb,coment write(*,'(1x,a)')'X,Y,Z /Angstr:' do 2 i=na+1,na+nb read(2,*)(x0(i,j),j=1,3) C C...iteration of J within loop 380 rather than direct use of index K c avoids optimisation out of this loop by Microsoft Powerstation Fortran BACKSPACE(2) READ(2,'(A)')LINE NDOTS=0 J=0 DO 380 K=1,80 J=J+1 IF(LINE(J:J).EQ.'.')NDOTS=NDOTS+1 IF(NDOTS.EQ.3)GOTO 381 380 CONTINUE 381 jj=j DO 382 k=J+1,80 jj=jj+1 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) c write(3,372)(x0(i,j),j=1,3),label(i) write(3,370)q0(i) write(3,370)(q1(i,j),j=1,3) write(3,370)(q2(i,j),j=1,5) write(3,370)(q3(i,j),j=1,7) write(3,370)(q4(i,j),j=1,9) write(*,372)(x0(i,j)*bohr,j=1,3),label(i) 2 continue c if(iangst.ne.1.and.iangst.ne.0)iangst=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 302 CONTINUE 1257 close(2) close(3) c c c...convert spherical to cartesian multipoles c 600 CALL CONVER c c...initial orientation of B from from keyboard c 352 write(*,506) * ' Rotation of B -> type theta1, phi1, alpha (in degrees)', * ' -> or 555,,, to terminate', * ' -> or 0,0,0 if only translation', * ' ... ' 506 format(1x/a/a/a//a,$) read(*,*,err=352)theta,psi,alpha c if(theta.eq.555.)then open(3,file=filnam,access='APPEND') write(*,'(1x/'' Final Cartesians /Angstr:''/)') write(3,'(1x/'' Final Cartesians /Angstr:''/)') do 1002 i=na+1,na+nb write(*,372)(x(i,j)*bohr,j=1,3),label(i) write(3,372)(x(i,j)*bohr,j=1,3),label(i) 1002 continue write(*,'(1x/'' Final Cartesians /bohr: ''/)') write(3,'(1x/'' Final Cartesians /bohr: ''/)') do 1003 i=na+1,na+nb write(*,372)(x(i,j) ,j=1,3),label(i) write(3,372)(x(i,j) ,j=1,3),label(i) 1003 continue write(*,'(1x/)') close(3) stop endif c c...conversion to bohr and of B vector to polar coordinates (if necessaary) c 254 ifit=0 r(1)=0.d0 r(2)=0.d0 r(3)=0.d0 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)(par(i)/conv,i=4,6) write(*,'(1x,a)')'X,Y,Z /Angstr:' open(3,file=filnam,access='APPEND') write(3,199) write(3,281)(par(i)/conv,i=4,6) close(3) 281 format(1x/' ROTATION OF B:'/ * /' Theta, Psi, Alpha: ',3f15.6/) c c c...ROTATION, BACKCONVERSION AND OUTPUT c c CALL ROTATB(par) CALL CONVSP c open(3,file=filnam,access='APPEND') WRITE(3,'(1X/A/)')' ROTATED DMA OF MOLECULE B:' write(3,221)nb,coment do 371 i=na+1,na+nb write(*,372)(x(i,j)*bohr,j=1,3),label(i) write(3,372)(x(i,j),j=1,3),label(i) write(3,370)q0(i) write(3,370)(q1(i,j),j=1,3) write(3,370)(q2(i,j),j=1,5) write(3,370)(q3(i,j),j=1,7) write(3,370)(q4(i,j),j=1,9) do 371 j=1,3 x0(i,j)=x(i,j) 371 continue 370 format(9f10.6) 372 format(3f10.6,a) write(3,199) 199 format(1x/) close(3) c c c...TRANSLATIONS c c 601 write(*,500) 500 FORMAT(/1x, * ' Translations 0 = none'/ * 1x,' 1 = simple shift of Cartesians'/ * 1x,' 2 = shift of multipole origin'/ * 1x,' 555 = terminate'// * 1x,' ... ',$) read(*,'(i5)')itrans IF(ITRANS.EQ.0)GOTO 600 c if(itrans.eq.555)then open(3,file=filnam,access='APPEND') write(*,'(1x/'' Final Cartesians /Angstr:''/)') write(3,'(1x/'' Final Cartesians /Angstr:''/)') do 1000 i=na+1,na+nb write(*,372)(x(i,j)*bohr,j=1,3),label(i) write(3,372)(x(i,j)*bohr,j=1,3),label(i) 1000 continue write(*,'(1x/'' Final Cartesians /bohr:''/)') do 1001 i=na+1,na+nb write(*,372)(x(i,j) ,j=1,3),label(i) write(3,372)(x(i,j) ,j=1,3),label(i) 1001 continue write(*,'(1x/)') close(3) stop endif c if(itrans.ne.1.and.itrans.ne.2)GOTO 601 write(*,'(/1x,a,$)')' DeltaX, DeltaY, DeltaZ /Angstr: ' read(*,*,err=601)delx,dely,delz delx=delx/bohr dely=dely/bohr delz=delz/bohr c c...simple shift of Cartesians c if(itrans.eq.1)then do 501 i=na+1,na+nb x(i,1)=x(i,1)+delx x(i,2)=x(i,2)+dely x(i,3)=x(i,3)+delz 501 continue goto 502 endif c c...translation of multipole origin c c...convert spherical to cartesian multipoles c CALL CONVER c c...shift the multipole of b c CALL SHIFT(-delx,-dely,-delz) c c...convert back to spherical multipoles c CALL CONVSP c c 502 open(3,file=filnam,access='APPEND') if(itrans.eq.1)then write(3,505) * ' SIMPLE TRANSLATION of Cartesians of B by (Angstr and bohr)', * bohr*delx,bohr*dely,bohr*delz,delx,dely,delz else write(3,505)' SHIFT OF MULTIPOLE origin of B by:', * bohr*delx,bohr*dely,bohr*delz,delx,dely,delz endif 505 format(1x/a,//' d.X, d.Y, d.Z: ',3f15.6/17x,3f15.6) c WRITE(3,'(1X/A/)')' TRANSLATED DMA OF MOLECULE B:' write(3,221)nb,coment do 481 i=na+1,na+nb write(3,372)(x(i,j),j=1,3),label(i) write(3,370)q0(i) write(3,370)(q1(i,j),j=1,3) write(3,370)(q2(i,j),j=1,5) write(3,370)(q3(i,j),j=1,7) write(3,370)(q4(i,j),j=1,9) 481 continue write(3,199) close(3) C GOTO 600 c stop end c c_____________________________________________________________________________ c SUBROUTINE CONVSP c c Conversion from cartesian to spherical 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) c one3=1.d0/3.d0 c do 52 i=na+1,na+nb q1(i,1)=d(i,3) q1(i,2)=d(i,1) q1(i,3)=d(i,2) 52 continue c do 51 i=na+1,na+nb q2(i,1)=t(i,3,3) q2(i,2)=dsqrt(4.d0/3.d0)*t(i,1,3) q2(i,3)=dsqrt(4.d0/3.d0)*t(i,2,3) q2(i,4)=dsqrt(1.d0/3.d0)*(t(i,1,1)-t(i,2,2)) q2(i,5)=dsqrt(4.d0/3.d0)*t(i,1,2) 51 continue c do 50 i=na+1,na+nb q3(i,1)=u(i,3,3,3) q3(i,2)=dsqrt(3.d0/2.d0)*u(i,1,3,3) q3(i,3)=dsqrt(3.d0/2.d0)*u(i,2,3,3) q3(i,4)=dsqrt(3.d0/5.d0)*(u(i,1,1,3)-u(i,2,2,3)) q3(i,5)=2.d0*dsqrt(3.d0/5.d0)*u(i,1,2,3) q3(i,6)=dsqrt(1.d0/10.d0)*(u(i,1,1,1)-3.d0*u(i,1,2,2)) q3(i,7)=dsqrt(1.d0/10.d0)*(3.d0*u(i,1,1,2)-u(i,2,2,2)) 50 continue c do 53 i=na+1,na+nb q4(i,1)=p(i,3,3,3,3) q4(i,2)=dsqrt(8.d0/5.d0)*p(i,1,3,3,3) q4(i,3)=dsqrt(8.d0/5.d0)*p(i,2,3,3,3) q4(i,4)=dsqrt(4.d0/5.d0)*(p(i,1,1,3,3)-p(i,2,2,3,3)) q4(i,5)=dsqrt(16.d0/5.d0)*p(i,1,2,3,3) q4(i,6)=dsqrt(8.d0/35.d0)*(p(i,1,1,1,3)-3.d0*p(i,1,2,2,3)) q4(i,7)=dsqrt(8.d0/35.d0)*(3.d0*p(i,1,1,2,3)-p(i,2,2,2,3)) q4(i,8)=dsqrt(1.d0/35.d0)*(p(i,1,1,1,1)-6.d0*p(i,1,1,2,2) * +p(i,2,2,2,2)) q4(i,9)=dsqrt(16.d0/35.d0)*(p(i,1,1,1,2)-p(i,1,2,2,2)) 53 continue 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 SHIFT(x,y,z) c c Shift the origin of the multipole from X,Y,Z (p.114, notebook 3) 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),xx(50,3),qq(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) integer*2 i c xsq=x*x ysq=y*y zsq=z*z rsq=xsq+ysq+zsq c32=3.0d0/2.0d0 c53=5.0d0/3.0d0 c83=8.0d0/3.0d0 c23=2.0d0/3.0d0 c52=5.0d0/2.0d0 do 1 i=na+1,na+nb q =qq(i) dx =d(i,1) dy =d(i,2) dz =d(i,3) txx =t(i,1,1) tyy =t(i,2,2) tzz =t(i,3,3) txy =t(i,1,2) txz =t(i,1,3) tyz =t(i,2,3) oxxx=u(i,1,1,1) oyyy=u(i,2,2,2) ozzz=u(i,3,3,3) oxxy=u(i,1,1,2) oxxz=u(i,1,1,3) oxyy=u(i,1,2,2) oyyz=u(i,2,2,3) oxzz=u(i,1,3,3) oyzz=u(i,2,3,3) oxyz=u(i,1,2,3) c pdx =dx+q*x pdy =dy+q*y pdz =dz+q*z ptxx =txx+2.d0*x*dx-y*dy-z*dz+0.5d0*q*(3.d0*xsq-rsq) ptyy =tyy-x*dx+2.d0*y*dy-z*dz+0.5d0*q*(3.d0*ysq-rsq) ptzz =tzz-x*dx-y*dy+2.d0*z*dz+0.5d0*q*(3.d0*zsq-rsq) ptxy =txy+c32*y*dx+c32*x*dy+c32*x*y*q ptxz =txz+c32*z*dx+c32*x*dz+c32*x*z*q ptyz =tyz+c32*z*dy+c32*y*dz+c32*y*z*q poxxx=oxxx+3.d0*x*txx-2.d0*y*txy-2.d0*z*txz * -c32*dx*(rsq-3.d0*xsq)-3.d0*x*y*dy-3.d0*x*z*dz * -0.5d0*q*x*(3.d0*rsq-5.d0*xsq) poyyy=oyyy-2.d0*x*txy+3.d0*y*tyy-2.d0*z*tyz-3.d0*x*y*dx * -c32*dy*(rsq-3.d0*ysq)-3.d0*y*z*dz * -0.5d0*q*y*(3.d0*rsq-5.d0*ysq) pozzz=ozzz-2.d0*x*txz-2.d0*y*tyz+3.d0*z*tzz-3.d0*x*z*dx * -3.d0*y*z*dy-c32*dz*(rsq-3.d0*zsq) * -0.5d0*q*z*(3.d0*rsq-5.d0*zsq) poxxy=oxxy+c53*y*txx+c83*x*txy-c23*y*tyy-c23*z*tyz+4.d0*x*y*dx * -0.5d0*dy*(rsq-5.d0*xsq+2.d0*ysq)-y*z*dz * -0.5d0*q*y*(rsq-5.d0*xsq) poxxz=oxxz+c53*z*txx+c83*x*txz-c23*y*tyz-c23*z*tzz+4.d0*x*z*dx * -y*z*dy-0.5d0*dz*(rsq-5.d0*xsq+2.d0*zsq) * -0.5d0*q*z*(rsq-5.d0*xsq) poxyy=oxyy-c23*x*txx+c83*y*txy-c23*z*txz+c53*x*tyy * -0.5d0*dx*(rsq+2.d0*xsq-5.d0*ysq)+4.d0*x*y*dy * -x*z*dz-0.5d0*q*x*(rsq-5.d0*ysq) poyyz=oyyz-c23*x*txz+c53*z*tyy+c83*y*tyz-c23*z*tzz-x*z*dx * +4.d0*y*z*dy-0.5d0*dz*(rsq-5.d0*ysq+2.d0*zsq) * -0.5d0*q*z*(rsq-5.d0*ysq) poxzz=oxzz-c23*x*txx-c23*y*txy+c83*z*txz+c53*x*tzz * -0.5d0*dx*(rsq+2.d0*xsq-5.d0*zsq)-x*y*dy * +4.d0*x*z*dz-0.5d0*q*x*(rsq-5.d0*zsq) poyzz=oyzz-c23*x*txy-c23*y*tyy+c83*z*tyz+c53*y*tzz-x*y*dx * -0.5d0*dy*(rsq+2.d0*ysq-5.d0*zsq)+4.d0*y*z*dz * -0.5d0*q*y*(rsq-5.d0*zsq) poxyz=oxyz+c53*z*txy+c53*y*txz+c53*x*tyz+c52*y*z*dx * +c52*x*z*dy+c52*x*y*dz+c52*x*y*z*q c d(i,1) =pdx d(i,2) =pdy d(i,3) =pdz t(i,1,1) =ptxx t(i,2,2) =ptyy t(i,3,3) =ptzz t(i,1,2) =ptxy t(i,1,3) =ptxz t(i,2,3) =ptyz u(i,1,1,1)=poxxx u(i,2,2,2)=poyyy u(i,3,3,3)=pozzz u(i,1,1,2)=poxxy u(i,1,1,3)=poxxz u(i,1,2,2)=poxyy u(i,2,2,3)=poyyz u(i,1,3,3)=poxzz u(i,2,3,3)=poyzz u(i,1,2,3)=poxyz c xx(i,1)=xx(i,1)-x xx(i,2)=xx(i,2)-y xx(i,3)=xx(i,3)-z c 1 continue c return end c 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 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_____________________________________________________________________________