diff --git a/source/analyze.f b/source/analyze.f index 43295c2c3..edc19bc99 100644 --- a/source/analyze.f +++ b/source/analyze.f @@ -1183,12 +1183,14 @@ subroutine paramyze (active) header = .false. write (iout,740) 740 format (/,' Dipole Polarizability Parameters :', - & //,10x,'Atom Number',5x,'Alpha',5x,'Damp', + & //,9x,'Atom Number',4x,'Alpha xx',2x, + & 'Alpha yy',2x,'Alpha zz',3x,'Damp', & 6x,'Polarization Group',/) end if - write (iout,750) i,ia,polarity(i),thole(i), + write (iout,750) i,ia,polarity(1,i),polarity(2,i), + & polarity(3,i), thole(i), & (ip11(j,ia),j=1,np11(ia)) - 750 format (i6,3x,i6,6x,f10.4,f9.3,3x,120i6) + 750 format (i6,3x,i6,6x,3f10.4,f9.3,3x,20i6) end if end do end if @@ -2403,10 +2405,12 @@ subroutine amberyze (active) header = .false. write (iout,730) 730 format (/,' Dipole Polarizability Parameters :', - & //,10x,'Atom Number',9x,'Alpha',8x, + & //,10x,'Atom Number',9x,'Alpha xx',8x, + & 'Alpha yy',4x, 'Alpha zz',4x, & 'Polarization Group',/) end if - write (iout,740) i,ia,polarity(i), + write (iout,740) i,ia,polarity(1,i),polarity(2,i), + & polarity(3,i), & (ip11(j,ia),j=1,np11(ia)) 740 format (i6,3x,i6,10x,f10.4,5x,20i6) end if diff --git a/source/epolar.f b/source/epolar.f index 2c2056cb8..6d6a65dfc 100644 --- a/source/epolar.f +++ b/source/epolar.f @@ -1628,7 +1628,7 @@ subroutine epolar0e use potent use units implicit none - integer i,j,ii + integer i,j,ii,m real*8 e,f,fi,term real*8 xd,yd,zd real*8 xu,yu,zu @@ -1660,20 +1660,17 @@ subroutine epolar0e c OpenMP directives for the major loop structure c !$OMP PARALLEL default(private) -!$OMP& shared(npole,polarity,f,uind,udirp,ep) +!$OMP& shared(npole,rpolarityinv,f,uind,udirp,ep) !$OMP DO reduction(+:ep) schedule(guided) c c get polarization energy via induced dipoles times field c do i = 1, npole - if (polarity(i) .ne. 0.0d0) then - fi = f / polarity(i) - e = 0.0d0 - do j = 1, 3 - e = fi * uind(j,i) * udirp(j,i) - ep = ep + e + do j = 1, 3 + do m = 1, 3 + ep = ep + f * uind(j,i) * udirp(m,i) * rpolarityinv(m,j,i) end do - end if + end do end do c c OpenMP directives for the major loop structure @@ -1759,7 +1756,7 @@ subroutine eprecip integer m1,m2,m3 integer ntot,nff integer nf1,nf2,nf3 - real*8 e,r1,r2,r3 + real*8 e,r1,r2,r3,f real*8 h1,h2,h3 real*8 volterm,denom real*8 hsq,expterm @@ -1769,6 +1766,7 @@ subroutine eprecip real*8, allocatable :: fuind(:,:) real*8, allocatable :: fuinp(:,:) c + f = electric / dielec c c return if the Ewald coefficient is zero c @@ -1890,7 +1888,7 @@ subroutine eprecip call fphi_mpole (fphi) do i = 1, npole do j = 1, 20 - fphi(j,i) = electric * fphi(j,i) + fphi(j,i) = f * fphi(j,i) end do end do end if @@ -1926,7 +1924,7 @@ subroutine eprecip if (.not. use_bounds) then expterm = 0.5d0 * pi / xbox struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 - e = 0.5d0 * electric * expterm * struc2 + e = 0.5d0 * f * expterm * struc2 ep = ep + e end if c @@ -1939,6 +1937,7 @@ subroutine eprecip e = e + fuind(k,i)*fphi(k+1,i) end do end do +c e = 0.5d0 * f * e e = 0.5d0 * e ep = ep + e c diff --git a/source/epolar1.f b/source/epolar1.f index bd187c08a..7cb5e27fd 100644 --- a/source/epolar1.f +++ b/source/epolar1.f @@ -72,6 +72,8 @@ subroutine epolar1a use virial implicit none integer i,j,k,m + integer ix,iy,iz + integer kx,ky,kz integer ii,kk,jcell integer iax,iay,iaz real*8 e,f,damp,expdamp @@ -124,6 +126,9 @@ subroutine epolar1a real*8 rc3(3),rc5(3),rc7(3) real*8 trq(3),fix(3) real*8 fiy(3),fiz(3) + real*8 uikx,uiky,uikz + real*8 ttmi(3),ttmk(3) + real*8, allocatable :: tem(:,:) real*8, allocatable :: pscale(:) real*8, allocatable :: dscale(:) real*8, allocatable :: uscale(:) @@ -161,6 +166,7 @@ subroutine epolar1a allocate (uscale(n)) allocate (ufld(3,n)) allocate (dufld(6,n)) + allocate (tem(3,n)) c c set exclusion coefficients and arrays to store fields c @@ -174,6 +180,9 @@ subroutine epolar1a do j = 1, 6 dufld(j,i) = 0.0d0 end do + do j = 1, 3 + tem(j,i) = 0.0d0 + end do end do c c set conversion factor, cutoff and switching coefficients @@ -558,6 +567,90 @@ subroutine epolar1a end do end if c +c add anisotropic dipole torques +c + + ! Note: this "if" check is not needed; the extra terms computed will + ! be zero in the isotropic case. Right now the go-to is isotropic dipoles + ! so I added it in, even though these few terms will not be a bottleneck + if(is_aniso(i) .or. is_aniso(k)) then + uikx = uky*uizp + ukyp*uiz - ukz*uiyp - ukzp*uiy + uiky = ukz*uixp + ukzp*uix - ukx*uizp - ukxp*uiz + uikz = ukx*uiyp + ukxp*uiy - uky*uixp - ukyp*uix + ! Site i + term1 = -ck*dsr3 + drk*dsr5 - qrrk*dsr7 + term2 = -ck*psr3 + drk*psr5 - qrrk*psr7 + term3 = 2.0d0 * dsr5 + term4 = 2.0d0 * psr5 + tem(1,i) = tem(1,i) + psr3*(dky*uiz - dkz*uiy) + tem(2,i) = tem(2,i) + psr3*(dkz*uix - dkx*uiz) + tem(3,i) = tem(3,i) + psr3*(dkx*uiy - dky*uix) + tem(1,i) = tem(1,i) + dsr3*(dky*uizp - dkz*uiyp) + tem(2,i) = tem(2,i) + dsr3*(dkz*uixp - dkx*uizp) + tem(3,i) = tem(3,i) + dsr3*(dkx*uiyp - dky*uixp) + tem(1,i) = tem(1,i) + term1*(uiyp*zr - uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr - uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr - uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr - uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr - uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr - uiy*xr) + tem(1,i) = tem(1,i) + term3*(uiyp*qrkz - uizp*qrky) + tem(2,i) = tem(2,i) + term3*(uizp*qrkx - uixp*qrkz) + tem(3,i) = tem(3,i) + term3*(uixp*qrky - uiyp*qrkx) + tem(1,i) = tem(1,i) + term4*(uiy*qrkz - uiz*qrky) + tem(2,i) = tem(2,i) + term4*(uiz*qrkx - uix*qrkz) + tem(3,i) = tem(3,i) + term4*(uix*qrky - uiy*qrkx) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*urk*sr5 + term2 = uscale(kk)*urkp*sr5 + tem(1,i) = tem(1,i) + uscale(kk)*sr3*uikx + tem(2,i) = tem(2,i) + uscale(kk)*sr3*uiky + tem(3,i) = tem(3,i) + uscale(kk)*sr3*uikz + tem(1,i) = tem(1,i) + term1*(uiyp*zr-uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr-uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr-uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr-uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr-uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr-uiy*xr) + endif + ! Site k + term1 = ci*dsr3 + dri*dsr5 + qrri*dsr7 + term2 = ci*psr3 + dri*psr5 + qrri*psr7 + term3 = -2.0d0 * dsr5 + term4 = -2.0d0 * psr5 + tem(1,k) = tem(1,k) + dsr3*(diy*ukzp - diz*ukyp) + tem(2,k) = tem(2,k) + dsr3*(diz*ukxp - dix*ukzp) + tem(3,k) = tem(3,k) + dsr3*(dix*ukyp - diy*ukxp) + tem(1,k) = tem(1,k) + psr3*(diy*ukz - diz*uky) + tem(2,k) = tem(2,k) + psr3*(diz*ukx - dix*ukz) + tem(3,k) = tem(3,k) + psr3*(dix*uky - diy*ukx) + tem(1,k) = tem(1,k) + term1*(ukyp*zr - ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr - ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr - ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr - ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr - ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr - uky*xr) + tem(1,k) = tem(1,k) + term3*(ukyp*qriz - ukzp*qriy) + tem(2,k) = tem(2,k) + term3*(ukzp*qrix - ukxp*qriz) + tem(3,k) = tem(3,k) + term3*(ukxp*qriy - ukyp*qrix) + tem(1,k) = tem(1,k) + term4*(uky*qriz - ukz*qriy) + tem(2,k) = tem(2,k) + term4*(ukz*qrix - ukx*qriz) + tem(3,k) = tem(3,k) + term4*(ukx*qriy - uky*qrix) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*uri*sr5 + term2 = uscale(kk)*urip*sr5 + tem(1,k) = tem(1,k) - uscale(kk)*sr3*uikx + tem(2,k) = tem(2,k) - uscale(kk)*sr3*uiky + tem(3,k) = tem(3,k) - uscale(kk)*sr3*uikz + tem(1,k) = tem(1,k) + term1*(ukyp*zr-ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr-ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr-ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr-ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr-ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr-uky*xr) + endif + endif +c c increment gradient and virial due to Cartesian forces c dep(1,ii) = dep(1,ii) + frcx @@ -1062,6 +1155,86 @@ subroutine epolar1a dsr7 = 0.5d0 * dsr7 end if c +c add anisotropic dipole torques +c + if(is_aniso(i) .or. is_aniso(k)) then + uikx = uky*uizp + ukyp*uiz - ukz*uiyp - ukzp*uiy + uiky = ukz*uixp + ukzp*uix - ukx*uizp - ukxp*uiz + uikz = ukx*uiyp + ukxp*uiy - uky*uixp - ukyp*uix + ! Site i + term1 = -ck*dsr3 + drk*dsr5 - qrrk*dsr7 + term2 = -ck*psr3 + drk*psr5 - qrrk*psr7 + term3 = 2.0d0 * dsr5 + term4 = 2.0d0 * psr5 + tem(1,i) = tem(1,i) + psr3*(dky*uiz - dkz*uiy) + tem(2,i) = tem(2,i) + psr3*(dkz*uix - dkx*uiz) + tem(3,i) = tem(3,i) + psr3*(dkx*uiy - dky*uix) + tem(1,i) = tem(1,i) + dsr3*(dky*uizp - dkz*uiyp) + tem(2,i) = tem(2,i) + dsr3*(dkz*uixp - dkx*uizp) + tem(3,i) = tem(3,i) + dsr3*(dkx*uiyp - dky*uixp) + tem(1,i) = tem(1,i) + term1*(uiyp*zr - uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr - uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr - uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr - uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr - uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr - uiy*xr) + tem(1,i) = tem(1,i) + term3*(uiyp*qrkz - uizp*qrky) + tem(2,i) = tem(2,i) + term3*(uizp*qrkx - uixp*qrkz) + tem(3,i) = tem(3,i) + term3*(uixp*qrky - uiyp*qrkx) + tem(1,i) = tem(1,i) + term4*(uiy*qrkz - uiz*qrky) + tem(2,i) = tem(2,i) + term4*(uiz*qrkx - uix*qrkz) + tem(3,i) = tem(3,i) + term4*(uix*qrky - uiy*qrkx) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*urk*sr5 + term2 = uscale(kk)*urkp*sr5 + tem(1,i) = tem(1,i) + uscale(kk)*sr3*uikx + tem(2,i) = tem(2,i) + uscale(kk)*sr3*uiky + tem(3,i) = tem(3,i) + uscale(kk)*sr3*uikz + tem(1,i) = tem(1,i) + term1*(uiyp*zr-uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr-uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr-uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr-uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr-uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr-uiy*xr) + endif + ! Site k + term1 = ci*dsr3 + dri*dsr5 + qrri*dsr7 + term2 = ci*psr3 + dri*psr5 + qrri*psr7 + term3 = -2.0d0 * dsr5 + term4 = -2.0d0 * psr5 + tem(1,k) = tem(1,k) + dsr3*(diy*ukzp - diz*ukyp) + tem(2,k) = tem(2,k) + dsr3*(diz*ukxp - dix*ukzp) + tem(3,k) = tem(3,k) + dsr3*(dix*ukyp - diy*ukxp) + tem(1,k) = tem(1,k) + psr3*(diy*ukz - diz*uky) + tem(2,k) = tem(2,k) + psr3*(diz*ukx - dix*ukz) + tem(3,k) = tem(3,k) + psr3*(dix*uky - diy*ukx) + tem(1,k) = tem(1,k) + term1*(ukyp*zr - ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr - ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr - ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr - ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr - ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr - uky*xr) + tem(1,k) = tem(1,k) + term3*(ukyp*qriz - ukzp*qriy) + tem(2,k) = tem(2,k) + term3*(ukzp*qrix - ukxp*qriz) + tem(3,k) = tem(3,k) + term3*(ukxp*qriy - ukyp*qrix) + tem(1,k) = tem(1,k) + term4*(uky*qriz - ukz*qriy) + tem(2,k) = tem(2,k) + term4*(ukz*qrix - ukx*qriz) + tem(3,k) = tem(3,k) + term4*(ukx*qriy - uky*qrix) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*uri*sr5 + term2 = uscale(kk)*urip*sr5 + tem(1,k) = tem(1,k) - uscale(kk)*sr3*uikx + tem(2,k) = tem(2,k) - uscale(kk)*sr3*uiky + tem(3,k) = tem(3,k) - uscale(kk)*sr3*uikz + tem(1,k) = tem(1,k) + term1*(ukyp*zr-ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr-ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr-ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr-ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr-ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr-uky*xr) + endif + endif +c c increment gradient and virial due to Cartesian forces c dep(1,ii) = dep(1,ii) + frcx @@ -1192,6 +1365,10 @@ subroutine epolar1a & + qiyz*dufld(4,i) - qixz*dufld(5,i) & + 2.0d0*qixy*(dufld(1,i)-dufld(3,i)) & + (qiyy-qixx)*dufld(2,i) + ! Anisotropic dipole terms + trq(1) = trq(1) + tem(1,i) + trq(2) = trq(2) + tem(2,i) + trq(3) = trq(3) + tem(3,i) call torque (i,trq,fix,fiy,fiz,dep) ii = ipole(i) iaz = zaxis(i) @@ -1233,6 +1410,7 @@ subroutine epolar1a deallocate (uscale) deallocate (ufld) deallocate (dufld) + deallocate (tem) return end c @@ -1321,6 +1499,9 @@ subroutine epolar1b real*8 rc3(3),rc5(3),rc7(3) real*8 trq(3),fix(3) real*8 fiy(3),fiz(3) + real*8 uikx,uiky,uikz + real*8 ttmi(3),ttmk(3) + real*8, allocatable :: tem(:,:) real*8, allocatable :: pscale(:) real*8, allocatable :: dscale(:) real*8, allocatable :: uscale(:) @@ -1358,6 +1539,7 @@ subroutine epolar1b allocate (uscale(n)) allocate (ufld(3,n)) allocate (dufld(6,n)) + allocate (tem(3,n)) c c set exclusion coefficients and arrays to store fields c @@ -1371,6 +1553,9 @@ subroutine epolar1b do j = 1, 6 dufld(j,i) = 0.0d0 end do + do j = 1, 3 + tem(j,i) = 0.0d0 + end do end do c c set conversion factor, cutoff and switching coefficients @@ -1382,6 +1567,7 @@ subroutine epolar1b c OpenMP directives for the major loop structure c !$OMP PARALLEL default(private) shared(npole,ipole,pdamp,thole, +!$OMP& is_aniso,tem, !$OMP& x,y,z,xaxis,yaxis,zaxis,rpole,uind,uinp,n12,i12,n13,i13,n14, !$OMP& i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14,p2scale, !$OMP& p3scale,p4scale,p41scale,p5scale,d1scale,d2scale,d3scale, @@ -1389,7 +1575,7 @@ subroutine epolar1b !$OMP& off2,f,molcule,coptmax,copm,uopt,uoptp,poltyp) !$OMP& shared (ep,einter,dep,vir,ufld,dufld) !$OMP& firstprivate(pscale,dscale,uscale) -!$OMP DO reduction(+:ep,einter,dep,vir,ufld,dufld) schedule(guided) +!$OMP DO reduction(+:ep,einter,dep,vir,ufld,dufld,tem) schedule(guided) c c compute the dipole polarization gradient components c @@ -1768,6 +1954,86 @@ subroutine epolar1b end do end if c +c add anisotropic dipole torques +c + if(is_aniso(i) .or. is_aniso(k)) then + uikx = uky*uizp + ukyp*uiz - ukz*uiyp - ukzp*uiy + uiky = ukz*uixp + ukzp*uix - ukx*uizp - ukxp*uiz + uikz = ukx*uiyp + ukxp*uiy - uky*uixp - ukyp*uix + ! Site i + term1 = -ck*dsr3 + drk*dsr5 - qrrk*dsr7 + term2 = -ck*psr3 + drk*psr5 - qrrk*psr7 + term3 = 2.0d0 * dsr5 + term4 = 2.0d0 * psr5 + tem(1,i) = tem(1,i) + psr3*(dky*uiz - dkz*uiy) + tem(2,i) = tem(2,i) + psr3*(dkz*uix - dkx*uiz) + tem(3,i) = tem(3,i) + psr3*(dkx*uiy - dky*uix) + tem(1,i) = tem(1,i) + dsr3*(dky*uizp - dkz*uiyp) + tem(2,i) = tem(2,i) + dsr3*(dkz*uixp - dkx*uizp) + tem(3,i) = tem(3,i) + dsr3*(dkx*uiyp - dky*uixp) + tem(1,i) = tem(1,i) + term1*(uiyp*zr - uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr - uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr - uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr - uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr - uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr - uiy*xr) + tem(1,i) = tem(1,i) + term3*(uiyp*qrkz - uizp*qrky) + tem(2,i) = tem(2,i) + term3*(uizp*qrkx - uixp*qrkz) + tem(3,i) = tem(3,i) + term3*(uixp*qrky - uiyp*qrkx) + tem(1,i) = tem(1,i) + term4*(uiy*qrkz - uiz*qrky) + tem(2,i) = tem(2,i) + term4*(uiz*qrkx - uix*qrkz) + tem(3,i) = tem(3,i) + term4*(uix*qrky - uiy*qrkx) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*urk*sr5 + term2 = uscale(kk)*urkp*sr5 + tem(1,i) = tem(1,i) + uscale(kk)*sr3*uikx + tem(2,i) = tem(2,i) + uscale(kk)*sr3*uiky + tem(3,i) = tem(3,i) + uscale(kk)*sr3*uikz + tem(1,i) = tem(1,i) + term1*(uiyp*zr-uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr-uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr-uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr-uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr-uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr-uiy*xr) + endif + ! Site k + term1 = ci*dsr3 + dri*dsr5 + qrri*dsr7 + term2 = ci*psr3 + dri*psr5 + qrri*psr7 + term3 = -2.0d0 * dsr5 + term4 = -2.0d0 * psr5 + tem(1,k) = tem(1,k) + dsr3*(diy*ukzp - diz*ukyp) + tem(2,k) = tem(2,k) + dsr3*(diz*ukxp - dix*ukzp) + tem(3,k) = tem(3,k) + dsr3*(dix*ukyp - diy*ukxp) + tem(1,k) = tem(1,k) + psr3*(diy*ukz - diz*uky) + tem(2,k) = tem(2,k) + psr3*(diz*ukx - dix*ukz) + tem(3,k) = tem(3,k) + psr3*(dix*uky - diy*ukx) + tem(1,k) = tem(1,k) + term1*(ukyp*zr - ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr - ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr - ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr - ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr - ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr - uky*xr) + tem(1,k) = tem(1,k) + term3*(ukyp*qriz - ukzp*qriy) + tem(2,k) = tem(2,k) + term3*(ukzp*qrix - ukxp*qriz) + tem(3,k) = tem(3,k) + term3*(ukxp*qriy - ukyp*qrix) + tem(1,k) = tem(1,k) + term4*(uky*qriz - ukz*qriy) + tem(2,k) = tem(2,k) + term4*(ukz*qrix - ukx*qriz) + tem(3,k) = tem(3,k) + term4*(ukx*qriy - uky*qrix) + if (poltyp .eq. 'MUTUAL') then + term1 = uscale(kk)*uri*sr5 + term2 = uscale(kk)*urip*sr5 + tem(1,k) = tem(1,k) - uscale(kk)*sr3*uikx + tem(2,k) = tem(2,k) - uscale(kk)*sr3*uiky + tem(3,k) = tem(3,k) - uscale(kk)*sr3*uikz + tem(1,k) = tem(1,k) + term1*(ukyp*zr-ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr-ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr-ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr-ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr-ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr-uky*xr) + endif + endif +c c increment gradient and virial due to Cartesian forces c dep(1,ii) = dep(1,ii) + frcx @@ -1901,6 +2167,10 @@ subroutine epolar1b & + qiyz*dufld(4,i) - qixz*dufld(5,i) & + 2.0d0*qixy*(dufld(1,i)-dufld(3,i)) & + (qiyy-qixx)*dufld(2,i) + ! Anisotropic dipole terms + trq(1) = trq(1) + tem(1,i) + trq(2) = trq(2) + tem(2,i) + trq(3) = trq(3) + tem(3,i) call torque (i,trq,fix,fiy,fiz,dep) ii = ipole(i) iaz = zaxis(i) @@ -1947,6 +2217,7 @@ subroutine epolar1b deallocate (uscale) deallocate (ufld) deallocate (dufld) + deallocate (tem) return end c @@ -2047,6 +2318,9 @@ subroutine epolar1c c term = (4.0d0/3.0d0) * f * aewald**3 / sqrtpi do i = 1, npole + ! anisotropic induced dipole torques are equal and opposite to + ! permanent dipole torques and, thus, they cancel each other + if(is_aniso(i)) cycle dix = rpole(2,i) diy = rpole(3,i) diz = rpole(4,i) @@ -2241,6 +2515,9 @@ subroutine epreal1c real*8 trq(3),fix(3) real*8 fiy(3),fiz(3) real*8 bn(0:4) + real*8 uikx,uiky,uikz + real*8 ttmi(3),ttmk(3) + real*8, allocatable :: tem(:,:) real*8, allocatable :: pscale(:) real*8, allocatable :: dscale(:) real*8, allocatable :: uscale(:) @@ -2257,6 +2534,7 @@ subroutine epreal1c allocate (uscale(n)) allocate (ufld(3,n)) allocate (dufld(6,n)) + allocate (tem(3,n)) c c set exclusion coefficients and arrays to store fields c @@ -2270,6 +2548,9 @@ subroutine epreal1c do j = 1, 6 dufld(j,i) = 0.0d0 end do + do j = 1, 3 + tem(j,i) = 0.0d0 + end do end do c c set conversion factor, cutoff and switching coefficients @@ -2761,6 +3042,86 @@ subroutine epreal1c end do end if c +c add anisotropic dipole torques +c + if(is_aniso(i) .or. is_aniso(k)) then + uikx = uky*uizp + ukyp*uiz - ukz*uiyp - ukzp*uiy + uiky = ukz*uixp + ukzp*uix - ukx*uizp - ukxp*uiz + uikz = ukx*uiyp + ukxp*uiy - uky*uixp - ukyp*uix + ! Site i + term1 = -ck*dsr3 + drk*dsr5 - qrrk*dsr7 + term2 = -ck*psr3 + drk*psr5 - qrrk*psr7 + term3 = 2.0d0 * dsr5 + term4 = 2.0d0 * psr5 + tem(1,i) = tem(1,i) + psr3*(dky*uiz - dkz*uiy) + tem(2,i) = tem(2,i) + psr3*(dkz*uix - dkx*uiz) + tem(3,i) = tem(3,i) + psr3*(dkx*uiy - dky*uix) + tem(1,i) = tem(1,i) + dsr3*(dky*uizp - dkz*uiyp) + tem(2,i) = tem(2,i) + dsr3*(dkz*uixp - dkx*uizp) + tem(3,i) = tem(3,i) + dsr3*(dkx*uiyp - dky*uixp) + tem(1,i) = tem(1,i) + term1*(uiyp*zr - uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr - uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr - uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr - uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr - uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr - uiy*xr) + tem(1,i) = tem(1,i) + term3*(uiyp*qrkz - uizp*qrky) + tem(2,i) = tem(2,i) + term3*(uizp*qrkx - uixp*qrkz) + tem(3,i) = tem(3,i) + term3*(uixp*qrky - uiyp*qrkx) + tem(1,i) = tem(1,i) + term4*(uiy*qrkz - uiz*qrky) + tem(2,i) = tem(2,i) + term4*(uiz*qrkx - uix*qrkz) + tem(3,i) = tem(3,i) + term4*(uix*qrky - uiy*qrkx) + if (poltyp .eq. 'MUTUAL') then + term1 = usr5*urk + term2 = usr5*urkp + tem(1,i) = tem(1,i) + usr3*uikx + tem(2,i) = tem(2,i) + usr3*uiky + tem(3,i) = tem(3,i) + usr3*uikz + tem(1,i) = tem(1,i) + term1*(uiyp*zr-uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr-uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr-uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr-uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr-uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr-uiy*xr) + endif + ! Site k + term1 = ci*dsr3 + dri*dsr5 + qrri*dsr7 + term2 = ci*psr3 + dri*psr5 + qrri*psr7 + term3 = -2.0d0 * dsr5 + term4 = -2.0d0 * psr5 + tem(1,k) = tem(1,k) + dsr3*(diy*ukzp - diz*ukyp) + tem(2,k) = tem(2,k) + dsr3*(diz*ukxp - dix*ukzp) + tem(3,k) = tem(3,k) + dsr3*(dix*ukyp - diy*ukxp) + tem(1,k) = tem(1,k) + psr3*(diy*ukz - diz*uky) + tem(2,k) = tem(2,k) + psr3*(diz*ukx - dix*ukz) + tem(3,k) = tem(3,k) + psr3*(dix*uky - diy*ukx) + tem(1,k) = tem(1,k) + term1*(ukyp*zr - ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr - ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr - ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr - ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr - ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr - uky*xr) + tem(1,k) = tem(1,k) + term3*(ukyp*qriz - ukzp*qriy) + tem(2,k) = tem(2,k) + term3*(ukzp*qrix - ukxp*qriz) + tem(3,k) = tem(3,k) + term3*(ukxp*qriy - ukyp*qrix) + tem(1,k) = tem(1,k) + term4*(uky*qriz - ukz*qriy) + tem(2,k) = tem(2,k) + term4*(ukz*qrix - ukx*qriz) + tem(3,k) = tem(3,k) + term4*(ukx*qriy - uky*qrix) + if (poltyp .eq. 'MUTUAL') then + term1 = usr5*uri + term2 = usr5*urip + tem(1,k) = tem(1,k) - usr3*uikx + tem(2,k) = tem(2,k) - usr3*uiky + tem(3,k) = tem(3,k) - usr3*uikz + tem(1,k) = tem(1,k) + term1*(ukyp*zr-ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr-ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr-ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr-ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr-ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr-uky*xr) + endif + endif +c c increment gradient and virial due to Cartesian forces c dep(1,ii) = dep(1,ii) - frcx @@ -3504,6 +3865,10 @@ subroutine epreal1c & + qiyz*dufld(4,i) - qixz*dufld(5,i) & + 2.0d0*qixy*(dufld(1,i)-dufld(3,i)) & + (qiyy-qixx)*dufld(2,i) + ! Anisotropic dipole terms + trq(1) = trq(1) + tem(1,i) + trq(2) = trq(2) + tem(2,i) + trq(3) = trq(3) + tem(3,i) call torque (i,trq,fix,fiy,fiz,dep) ii = ipole(i) iaz = zaxis(i) @@ -3545,6 +3910,7 @@ subroutine epreal1c deallocate (uscale) deallocate (ufld) deallocate (dufld) + deallocate (tem) return end c @@ -3645,6 +4011,9 @@ subroutine epolar1d c term = (4.0d0/3.0d0) * f * aewald**3 / sqrtpi do i = 1, npole + ! anisotropic induced dipole torques are equal and opposite to + ! permanent dipole torques and, thus, they cancel each other + if(is_aniso(i)) cycle dix = rpole(2,i) diy = rpole(3,i) diz = rpole(4,i) @@ -3839,6 +4208,9 @@ subroutine epreal1d real*8 trq(3),fix(3) real*8 fiy(3),fiz(3) real*8 bn(0:4) + real*8 uikx,uiky,uikz + real*8 ttmi(3),ttmk(3) + real*8, allocatable :: tem(:,:) real*8, allocatable :: pscale(:) real*8, allocatable :: dscale(:) real*8, allocatable :: uscale(:) @@ -3855,6 +4227,7 @@ subroutine epreal1d allocate (uscale(n)) allocate (ufld(3,n)) allocate (dufld(6,n)) + allocate (tem(3,n)) c c set exclusion coefficients and arrays to store fields c @@ -3868,6 +4241,9 @@ subroutine epreal1d do j = 1, 6 dufld(j,i) = 0.0d0 end do + do j = 1, 3 + tem(j,i) = 0.0d0 + end do end do c c set conversion factor, cutoff and switching coefficients @@ -3879,6 +4255,7 @@ subroutine epreal1d c OpenMP directives for the major loop structure c !$OMP PARALLEL default(private) shared(npole,ipole,pdamp,thole, +!$OMP& is_aniso,tem, !$OMP& x,y,z,xaxis,yaxis,zaxis,rpole,uind,uinp,n12,i12,n13,i13,n14, !$OMP& i14,n15,i15,np11,ip11,np12,ip12,np13,ip13,np14,ip14,p2scale, !$OMP& p3scale,p4scale,p41scale,p5scale,d1scale,d2scale,d3scale, @@ -3886,7 +4263,7 @@ subroutine epreal1d !$OMP& off2,f,aewald,molcule,coptmax,copm,uopt,uoptp,poltyp) !$OMP& shared (ep,einter,dep,vir,ufld,dufld) !$OMP& firstprivate(pscale,dscale,uscale) -!$OMP DO reduction(+:ep,einter,dep,vir,ufld,dufld) schedule(guided) +!$OMP DO reduction(+:ep,einter,dep,vir,ufld,dufld,tem) schedule(guided) c c compute the dipole polarization gradient components c @@ -4372,6 +4749,86 @@ subroutine epreal1d end do end if c +c add anisotropic dipole torques +c + if(is_aniso(i) .or. is_aniso(k)) then + uikx = uky*uizp + ukyp*uiz - ukz*uiyp - ukzp*uiy + uiky = ukz*uixp + ukzp*uix - ukx*uizp - ukxp*uiz + uikz = ukx*uiyp + ukxp*uiy - uky*uixp - ukyp*uix + ! Site i + term1 = -ck*dsr3 + drk*dsr5 - qrrk*dsr7 + term2 = -ck*psr3 + drk*psr5 - qrrk*psr7 + term3 = 2.0d0 * dsr5 + term4 = 2.0d0 * psr5 + tem(1,i) = tem(1,i) + psr3*(dky*uiz - dkz*uiy) + tem(2,i) = tem(2,i) + psr3*(dkz*uix - dkx*uiz) + tem(3,i) = tem(3,i) + psr3*(dkx*uiy - dky*uix) + tem(1,i) = tem(1,i) + dsr3*(dky*uizp - dkz*uiyp) + tem(2,i) = tem(2,i) + dsr3*(dkz*uixp - dkx*uizp) + tem(3,i) = tem(3,i) + dsr3*(dkx*uiyp - dky*uixp) + tem(1,i) = tem(1,i) + term1*(uiyp*zr - uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr - uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr - uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr - uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr - uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr - uiy*xr) + tem(1,i) = tem(1,i) + term3*(uiyp*qrkz - uizp*qrky) + tem(2,i) = tem(2,i) + term3*(uizp*qrkx - uixp*qrkz) + tem(3,i) = tem(3,i) + term3*(uixp*qrky - uiyp*qrkx) + tem(1,i) = tem(1,i) + term4*(uiy*qrkz - uiz*qrky) + tem(2,i) = tem(2,i) + term4*(uiz*qrkx - uix*qrkz) + tem(3,i) = tem(3,i) + term4*(uix*qrky - uiy*qrkx) + if (poltyp .eq. 'MUTUAL') then + term1 = usr5*urk + term2 = usr5*urkp + tem(1,i) = tem(1,i) + usr3*uikx + tem(2,i) = tem(2,i) + usr3*uiky + tem(3,i) = tem(3,i) + usr3*uikz + tem(1,i) = tem(1,i) + term1*(uiyp*zr-uizp*yr) + tem(2,i) = tem(2,i) + term1*(uizp*xr-uixp*zr) + tem(3,i) = tem(3,i) + term1*(uixp*yr-uiyp*xr) + tem(1,i) = tem(1,i) + term2*(uiy*zr-uiz*yr) + tem(2,i) = tem(2,i) + term2*(uiz*xr-uix*zr) + tem(3,i) = tem(3,i) + term2*(uix*yr-uiy*xr) + endif + ! Site k + term1 = ci*dsr3 + dri*dsr5 + qrri*dsr7 + term2 = ci*psr3 + dri*psr5 + qrri*psr7 + term3 = -2.0d0 * dsr5 + term4 = -2.0d0 * psr5 + tem(1,k) = tem(1,k) + dsr3*(diy*ukzp - diz*ukyp) + tem(2,k) = tem(2,k) + dsr3*(diz*ukxp - dix*ukzp) + tem(3,k) = tem(3,k) + dsr3*(dix*ukyp - diy*ukxp) + tem(1,k) = tem(1,k) + psr3*(diy*ukz - diz*uky) + tem(2,k) = tem(2,k) + psr3*(diz*ukx - dix*ukz) + tem(3,k) = tem(3,k) + psr3*(dix*uky - diy*ukx) + tem(1,k) = tem(1,k) + term1*(ukyp*zr - ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr - ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr - ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr - ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr - ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr - uky*xr) + tem(1,k) = tem(1,k) + term3*(ukyp*qriz - ukzp*qriy) + tem(2,k) = tem(2,k) + term3*(ukzp*qrix - ukxp*qriz) + tem(3,k) = tem(3,k) + term3*(ukxp*qriy - ukyp*qrix) + tem(1,k) = tem(1,k) + term4*(uky*qriz - ukz*qriy) + tem(2,k) = tem(2,k) + term4*(ukz*qrix - ukx*qriz) + tem(3,k) = tem(3,k) + term4*(ukx*qriy - uky*qrix) + if (poltyp .eq. 'MUTUAL') then + term1 = usr5*uri + term2 = usr5*urip + tem(1,k) = tem(1,k) - usr3*uikx + tem(2,k) = tem(2,k) - usr3*uiky + tem(3,k) = tem(3,k) - usr3*uikz + tem(1,k) = tem(1,k) + term1*(ukyp*zr-ukzp*yr) + tem(2,k) = tem(2,k) + term1*(ukzp*xr-ukxp*zr) + tem(3,k) = tem(3,k) + term1*(ukxp*yr-ukyp*xr) + tem(1,k) = tem(1,k) + term2*(uky*zr-ukz*yr) + tem(2,k) = tem(2,k) + term2*(ukz*xr-ukx*zr) + tem(3,k) = tem(3,k) + term2*(ukx*yr-uky*xr) + endif + endif +c c increment gradient and virial due to Cartesian forces c dep(1,ii) = dep(1,ii) - frcx @@ -4505,6 +4962,12 @@ subroutine epreal1d & + qiyz*dufld(4,i) - qixz*dufld(5,i) & + 2.0d0*qixy*(dufld(1,i)-dufld(3,i)) & + (qiyy-qixx)*dufld(2,i) +c +c anisotropic dipole terms +c + trq(1) = trq(1) + tem(1,i) + trq(2) = trq(2) + tem(2,i) + trq(3) = trq(3) + tem(3,i) call torque (i,trq,fix,fiy,fiz,dep) ii = ipole(i) iaz = zaxis(i) @@ -4606,6 +5069,7 @@ subroutine eprecip1 integer deriv2(10) integer deriv3(10) real*8 e,eterm + real*8 f real*8 r1,r2,r3 real*8 h1,h2,h3 real*8 f1,f2,f3 @@ -4615,6 +5079,7 @@ subroutine eprecip1 real*8 hsq,expterm real*8 term,pterm real*8 vterm,struc2 + real*8 idip(3) real*8 trq(3),fix(3) real*8 fiy(3),fiz(3) real*8 cphim(4),cphid(4) @@ -4638,6 +5103,10 @@ subroutine eprecip1 c if (aewald .lt. 1.0d-6) return c +c add dielectric +c + f = electric / dielec +c c perform dynamic allocation of some global arrays c if (allocated(cmp)) then @@ -4744,7 +5213,7 @@ subroutine eprecip1 if (mod(m1+m2+m3,2) .ne. 0) expterm = 0.0d0 end if struc2 = qgrid(1,k1,k2,k3)**2 + qgrid(2,k1,k2,k3)**2 - eterm = 0.5d0 * electric * expterm * struc2 + eterm = 0.5d0 * f * expterm * struc2 vterm = (2.0d0/hsq) * (1.0d0-term) * eterm vxx = vxx - h1*h1*vterm + eterm vxy = vxy - h1*h2*vterm @@ -4761,7 +5230,7 @@ subroutine eprecip1 if (.not. use_bounds) then expterm = 0.5d0 * pi / xbox struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 - e = 0.5d0 * expterm * struc2 + e = 0.5d0 * f * expterm * struc2 em = em + e qfac(1,1,1) = expterm end if @@ -4784,12 +5253,26 @@ subroutine eprecip1 call fphi_mpole (fphi) do i = 1, npole do j = 1, 20 - fphi(j,i) = electric * fphi(j,i) + fphi(j,i) = f * fphi(j,i) end do end do call fphi_to_cphi (fphi,cphi) end if c +c handle anisotropic induced dipoles +c + do i = 1, npole + if(is_aniso(i)) then + idip(1) = 0.5d0*(uind(1,i) + uinp(1,i)) + idip(2) = 0.5d0*(uind(2,i) + uinp(2,i)) + idip(3) = 0.5d0*(uind(3,i) + uinp(3,i)) + trq(1) = idip(3)*cphi(3,i) - idip(2)*cphi(4,i) + trq(2) = idip(1)*cphi(4,i) - idip(3)*cphi(2,i) + trq(3) = idip(2)*cphi(2,i) - idip(1)*cphi(3,i) + call torque (i,trq,fix,fiy,fiz,dep) + endif + end do +c c perform dynamic allocation of some local arrays c allocate (fuind(3,npole)) @@ -4824,7 +5307,7 @@ subroutine eprecip1 if (.not. use_bounds) then expterm = 0.5d0 * pi / xbox struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 - e = 0.5d0 * expterm * struc2 + e = 0.5d0 * f * expterm * struc2 ep = ep + e end if c @@ -4846,11 +5329,11 @@ subroutine eprecip1 call fphi_uind (fphid,fphip,fphidp) do i = 1, npole do j = 1, 10 - fphid(j,i) = electric * fphid(j,i) - fphip(j,i) = electric * fphip(j,i) + fphid(j,i) = f * fphid(j,i) + fphip(j,i) = f * fphip(j,i) end do do j = 1, 20 - fphidp(j,i) = electric * fphidp(j,i) + fphidp(j,i) = f * fphidp(j,i) end do end do c @@ -4927,6 +5410,15 @@ subroutine eprecip1 & + 2.0d0*(cmp(6,i)-cmp(5,i))*cphi(8,i) & + cmp(8,i)*cphi(5,i) + cmp(10,i)*cphi(9,i) & - cmp(8,i)*cphi(6,i) - cmp(9,i)*cphi(10,i) + if(is_aniso(i) .and. poltyp .eq. "MUTUAL") then + ! add mutual induced torques + idip(1) = 0.5d0*(uind(1,i) + uinp(1,i)) + idip(2) = 0.5d0*(uind(2,i) + uinp(2,i)) + idip(3) = 0.5d0*(uind(3,i) + uinp(3,i)) + trq(1) = trq(1) + idip(3)*cphi(3,i) - idip(2)*cphi(4,i) + trq(2) = trq(2) + idip(1)*cphi(4,i) - idip(3)*cphi(2,i) + trq(3) = trq(3) + idip(2)*cphi(2,i) - idip(1)*cphi(3,i) + endif call torque (i,trq,fix,fiy,fiz,dep) end do c @@ -5002,8 +5494,8 @@ subroutine eprecip1 ii = ipole(i) do k = 0, coptmax-1 do j = 1, 10 - fphid(j,i) = electric * fopt(k,j,i) - fphip(j,i) = electric * foptp(k,j,i) + fphid(j,i) = f * fopt(k,j,i) + fphip(j,i) = f * foptp(k,j,i) end do do m = 0, coptmax-k-1 do j = 1, 3 @@ -5140,7 +5632,7 @@ subroutine eprecip1 end if struc2 = qgrid(1,k1,k2,k3)*qgrip(1,k1,k2,k3) & + qgrid(2,k1,k2,k3)*qgrip(2,k1,k2,k3) - eterm = 0.5d0 * electric * expterm * struc2 + eterm = 0.5d0 * f * expterm * struc2 vterm = (2.0d0/hsq) * (1.0d0-term) * eterm vxx = vxx + h1*h1*vterm - eterm vxy = vxy + h1*h2*vterm @@ -5216,7 +5708,7 @@ subroutine eprecip1 end if struc2 = qgrid(1,k1,k2,k3)*qgrip(1,k1,k2,k3) & + qgrid(2,k1,k2,k3)*qgrip(2,k1,k2,k3) - eterm = 0.5d0 * electric * expterm * struc2 + eterm = 0.5d0 * f * expterm * struc2 vterm = (2.0d0/hsq) * (1.0d0-term) * eterm vxx = vxx - h1*h1*vterm + eterm vxy = vxy - h1*h2*vterm diff --git a/source/epolar3.f b/source/epolar3.f index 1d3219f59..f74654207 100644 --- a/source/epolar3.f +++ b/source/epolar3.f @@ -1923,7 +1923,7 @@ subroutine epolar3e use potent use units implicit none - integer i,j,ii + integer i,j,ii,m real*8 e,f,fi,term real*8 xd,yd,zd real*8 xu,yu,zu @@ -1977,12 +1977,15 @@ subroutine epolar3e c get polarization energy via induced dipoles times field c do i = 1, npole - if (polarity(i) .ne. 0.0d0) then - ii = ipole(i) - fi = f / polarity(i) + ii = ipole(i) + if (polarity(1,i) .ne. 0.0d0 .or. + & polarity(2,i) .ne. 0.0d0 .or. + & polarity(3,i) .ne. 0.0d0) then e = 0.0d0 do j = 1, 3 - e = e + fi*uind(j,i)*udirp(j,i) + do m = 1, 3 + e = e + f * uind(j,i) * udirp(m,i)*rpolarityinv(m,j,i) + end do end do nep = nep + 1 ep = ep + e @@ -1997,11 +2000,13 @@ subroutine epolar3e write (iout,20) 20 format (/,' Individual Dipole Polarization', & ' Interactions :', - & //,' Type',9x,'Atom Name',24x,'Alpha', + & //,' Type',9x,'Atom Name',24x,'Alpha xx', + & 'Alpha yy', 8x, 'Alpha zz', & 8x,'Energy',/) end if - write (iout,30) ii,name(ii),polarity(i),e - 30 format (' Polar',5x,i7,'-',a3,16x,2f14.4) + write (iout,30) ii,name(ii),polarity(1,i),polarity(2,i), + & polarity(3,i), e + 30 format (' Polar',5x,i7,'-',a3,16x,4f14.4) end if end if end do diff --git a/source/final.f b/source/final.f index 5b8a7cca2..6f60297d7 100644 --- a/source/final.f +++ b/source/final.f @@ -629,7 +629,10 @@ subroutine final c if (allocated(copt)) deallocate (copt) if (allocated(copm)) deallocate (copm) + if (allocated(is_aniso)) deallocate (is_aniso) if (allocated(polarity)) deallocate (polarity) + if (allocated(rpolarity)) deallocate (rpolarity) + if (allocated(rpolarityinv)) deallocate (rpolarityinv) if (allocated(thole)) deallocate (thole) if (allocated(pdamp)) deallocate (pdamp) if (allocated(udir)) deallocate (udir) diff --git a/source/induce.f b/source/induce.f index 00fb9e512..37db41f3b 100644 --- a/source/induce.f +++ b/source/induce.f @@ -76,7 +76,9 @@ subroutine induce if (debug .and. use_polar) then header = .true. do i = 1, npole - if (polarity(i) .ne. 0.0d0) then + if (polarity(1,i) .ne. 0.0d0 .or. + & polarity(2,i) .ne. 0.0d0 .or. + & polarity(3,i) .ne. 0.0d0) then if (header) then header = .false. if (solvtyp.eq.'GK' .or. solvtyp.eq.'PB') then @@ -104,13 +106,13 @@ subroutine induce k = ipole(i) norm = sqrt(uind(1,i)**2+uind(2,i)**2+uind(3,i)**2) if (digits .ge. 8) then - write (iout,60) k,(debye*uind(j,i),j=1,3),debye*norm + write (iout,60) k,(uind(j,i),j=1,3),debye*norm 60 format (i8,3x,4f16.8) else if (digits .ge. 6) then - write (iout,70) k,(debye*uind(j,i),j=1,3),debye*norm + write (iout,70) k,(uind(j,i),j=1,3),debye*norm 70 format (i8,4x,4f14.6) else - write (iout,80) k,(debye*uind(j,i),j=1,3),debye*norm + write (iout,80) k,(uind(j,i),j=1,3),debye*norm 80 format (i8,5x,4f12.4) end if end if @@ -118,7 +120,9 @@ subroutine induce header = .true. if (solvtyp.eq.'GK' .or. solvtyp.eq.'PB') then do i = 1, npole - if (polarity(i) .ne. 0.0d0) then + if (polarity(1,i) .ne. 0.0d0 .or. + & polarity(2,i) .ne. 0.0d0 .or. + & polarity(3,i) .ne. 0.0d0) then if (header) then header = .false. write (iout,90) @@ -186,15 +190,13 @@ subroutine induce0a use units use uprior implicit none - integer i,j,k,iter + integer i,j,k,m,iter integer maxiter - real*8 polmin real*8 eps,epsold real*8 epsd,epsp real*8 udsum,upsum real*8 a,ap,b,bp real*8 sum,sump - real*8, allocatable :: poli(:) real*8, allocatable :: field(:,:) real*8, allocatable :: fieldp(:,:) real*8, allocatable :: rsd(:,:) @@ -238,11 +240,16 @@ subroutine induce0a c c set induced dipoles to polarizability times direct field c + do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - udir(j,i) = polarity(i) * field(j,i) - udirp(j,i) = polarity(i) * fieldp(j,i) + udir(j,i) = 0.0d0 + udirp(j,i) = 0.0d0 + do m = 1, 3 + udir(j,i) = udir(j,i) + rpolarity(m,j,i)*field(m,i) + udirp(j,i) = udirp(j,i) + rpolarity(m,j,i)*fieldp(m,i) + end do uind(j,i) = udir(j,i) uinp(j,i) = udirp(j,i) end do @@ -272,8 +279,14 @@ subroutine induce0a do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uopt(k,j,i) = polarity(i) * field(j,i) - uoptp(k,j,i) = polarity(i) * fieldp(j,i) + uopt(k,j,i) = 0.0d0 + uoptp(k,j,i) = 0.0d0 + do m = 1, 3 + uopt(k,j,i) = uopt(k,j,i) + + & rpolarity(m,j,i) * field(m,i) + uoptp(k,j,i) = uoptp(k,j,i) + + & rpolarity(m,j,i) * fieldp(m,i) + end do uind(j,i) = uopt(k,j,i) uinp(j,i) = uoptp(k,j,i) end do @@ -308,7 +321,6 @@ subroutine induce0a done = .false. maxiter = 100 iter = 0 - polmin = 0.00000001d0 eps = 100.0d0 c c estimate induced dipoles using a polynomial predictor @@ -342,7 +354,6 @@ subroutine induce0a c c perform dynamic allocation of some local arrays c - allocate (poli(npole)) allocate (rsd(3,npole)) allocate (rsdp(3,npole)) allocate (zrsd(3,npole)) @@ -366,12 +377,15 @@ subroutine induce0a c do i = 1, npole if (douind(ipole(i))) then - poli(i) = max(polmin,polarity(i)) do j = 1, 3 - rsd(j,i) = (udir(j,i)-uind(j,i))/poli(i) - & + field(j,i) - rsdp(j,i) = (udirp(j,i)-uinp(j,i))/poli(i) - & + fieldp(j,i) + rsd(j,i) = field(j,i) + rsdp(j,i) = fieldp(j,i) + do m = 1, 3 + rsd(j,i) = rsd(j,i) + (udir(m,i)-uind(m,i)) + & * rpolarityinv(m,j,i) + rsdp(j,i) = rsdp(j,i) + (udirp(m,i)-uinp(m,i)) + & * rpolarityinv(m,j,i) + enddo end do else do j = 1, 3 @@ -425,8 +439,14 @@ subroutine induce0a do j = 1, 3 uind(j,i) = vec(j,i) uinp(j,i) = vecp(j,i) - vec(j,i) = conj(j,i)/poli(i) - field(j,i) - vecp(j,i) = conjp(j,i)/poli(i) - fieldp(j,i) + vec(j,i) = -field(j,i) + vecp(j,i) = -fieldp(j,i) + do m = 1, 3 + vec(j,i) = vec(j,i) + + & conj(m,i)*rpolarityinv(m,j,i) + vecp(j,i) = vecp(j,i) + + & conjp(m,i)*rpolarityinv(m,j,i) + end do end do end if end do @@ -511,8 +531,12 @@ subroutine induce0a do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uind(j,i) = uind(j,i) + poli(i)*rsd(j,i) - uinp(j,i) = uinp(j,i) + poli(i)*rsdp(j,i) + do m = 1, 3 + uind(j,i) = uind(j,i) + + & rpolarity(m,j,i)*rsd(m,i) + uinp(j,i) = uinp(j,i) + + & rpolarity(m,j,i)*rsdp(m,i) + end do end do end if end do @@ -521,7 +545,6 @@ subroutine induce0a c c perform deallocation of some local arrays c - deallocate (poli) deallocate (rsd) deallocate (rsdp) deallocate (zrsd) @@ -3433,9 +3456,8 @@ subroutine induce0d use units use uprior implicit none - integer i,j,k,iter + integer i,j,k,m,iter integer maxiter - real*8 polmin real*8 eps,epsold real*8 epsd,epsp real*8 epsds,epsps @@ -3445,7 +3467,6 @@ subroutine induce0d real*8 b,bp,bs,bps real*8 sum,sump real*8 sums,sumps - real*8, allocatable :: poli(:) real*8, allocatable :: field(:,:) real*8, allocatable :: fieldp(:,:) real*8, allocatable :: fields(:,:) @@ -3511,10 +3532,20 @@ subroutine induce0d do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - udir(j,i) = polarity(i) * field(j,i) - udirp(j,i) = polarity(i) * fieldp(j,i) - udirs(j,i) = polarity(i) * fields(j,i) - udirps(j,i) = polarity(i) * fieldps(j,i) + udir(j,i) = 0.0d0 + udirp(j,i) = 0.0d0 + udirs(j,i) = 0.0d0 + udirps(j,i) = 0.0d0 + do m = 1, 3 + udir(j,i) = udir(j,i) + rpolarity(m,j,i) + & * field(m,i) + udirp(j,i) = udirp(j,i) + rpolarity(m,j,i) + & * fieldp(m,i) + udirs(j,i) = udirs(j,i) + rpolarity(m,j,i) + & * fields(m,i) + udirps(j,i) = udirps(j,i) + rpolarity(m,j,i) + & * fieldps(m,i) + end do uind(j,i) = udir(j,i) uinp(j,i) = udirp(j,i) uinds(j,i) = udirs(j,i) @@ -3541,10 +3572,20 @@ subroutine induce0d do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uopt(k,j,i) = polarity(i) * field(j,i) - uoptp(k,j,i) = polarity(i) * fieldp(j,i) - uopts(k,j,i) = polarity(i) * fields(j,i) - uoptps(k,j,i) = polarity(i) * fieldps(j,i) + uopt(k,j,i) = 0.0d0 + uoptp(k,j,i) = 0.0d0 + uopts(k,j,i) = 0.0d0 + uoptps(k,j,i) = 0.0d0 + do m = 1, 3 + uopt(k,j,i) = uopt(k,j,i) + + & rpolarity(m,j,i) * field(m,i) + uoptp(k,j,i) = uoptp(k,j,i) + + & rpolarity(m,j,i) * fieldp(m,i) + uopts(k,j,i) = uopts(k,j,i) + + & rpolarity(m,j,i) * fields(m,i) + uoptps(k,j,i) = uoptps(k,j,i) + + & rpolarity(m,j,i) * fieldps(m,i) + end do uind(j,i) = uopt(k,j,i) uinp(j,i) = uoptp(k,j,i) uinds(j,i) = uopts(k,j,i) @@ -3593,7 +3634,6 @@ subroutine induce0d done = .false. maxiter = 100 iter = 0 - polmin = 0.00000001d0 eps = 100.0d0 c c estimated induced dipoles from polynomial predictor @@ -3621,7 +3661,6 @@ subroutine induce0d c c perform dynamic allocation of some local arrays c - allocate (poli(npole)) allocate (rsd(3,npole)) allocate (rsdp(3,npole)) allocate (rsds(3,npole)) @@ -3644,20 +3683,33 @@ subroutine induce0d call ufield0d (field,fieldp,fields,fieldps) do i = 1, npole if (douind(ipole(i))) then - poli(i) = max(polmin,polarity(i)) do j = 1, 3 - rsd(j,i) = (udir(j,i)-uind(j,i))/poli(i) - & + field(j,i) - rsdp(j,i) = (udirp(j,i)-uinp(j,i))/poli(i) - & + fieldp(j,i) - rsds(j,i) = (udirs(j,i)-uinds(j,i))/poli(i) - & + fields(j,i) - rsdps(j,i) = (udirps(j,i)-uinps(j,i))/poli(i) - & + fieldps(j,i) - zrsd(j,i) = rsd(j,i) * poli(i) - zrsdp(j,i) = rsdp(j,i) * poli(i) - zrsds(j,i) = rsds(j,i) * poli(i) - zrsdps(j,i) = rsdps(j,i) * poli(i) + rsd(j,i) = field(j,i) + rsdp(j,i) = fieldp(j,i) + rsds(j,i) = fields(j,i) + rsdps(j,i) = fieldps(j,i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + zrsds(j,i) = 0.0d0 + zrsdps(j,i) = 0.0d0 + do m = 1, 3 + rsd(j,i) = rsd(j,i) + + & (udir(m,i)-uind(m,i))*rpolarityinv(m,j,i) + rsdp(j,i) = rsdp(j,i) + + & (udirp(m,i)-uinp(m,i))*rpolarityinv(m,j,i) + rsds(j,i) = rsds(j,i) + + & (udirs(m,i)-uinds(m,i))*rpolarityinv(m,j,i) + rsdps(j,i) = rsdps(j,i) + + & (udirps(m,i)-uinps(m,i))*rpolarityinv(m,j,i) + zrsd(j,i) = zrsd(j,i) + + & rsd(m,i)*rpolarity(m,j,i) + zrsdp(j,i) = zrsdp(j,i) + + & rsdp(m,i)*rpolarity(m,j,i) + zrsds(j,i) = zrsds(j,i) + + & rsds(m,i)*rpolarity(m,j,i) + zrsdps(j,i) = zrsdps(j,i) + + & rsdps(m,i)*rpolarity(m,j,i) + end do conj(j,i) = zrsd(j,i) conjp(j,i) = zrsdp(j,i) conjs(j,i) = zrsds(j,i) @@ -3692,10 +3744,20 @@ subroutine induce0d uinp(j,i) = vecp(j,i) uinds(j,i) = vecs(j,i) uinps(j,i) = vecps(j,i) - vec(j,i) = conj(j,i)/poli(i) - field(j,i) - vecp(j,i) = conjp(j,i)/poli(i) - fieldp(j,i) - vecs(j,i) = conjs(j,i)/poli(i) - fields(j,i) - vecps(j,i) = conjps(j,i)/poli(i) - fieldps(j,i) + vec(j,i) = -field(j,i) + vecp(j,i) = -fieldp(j,i) + vecs(j,i) = -fields(j,i) + vecps(j,i) = -fieldps(j,i) + do m = 1, 3 + vec(j,i) = vec(j,i) + + & conj(m,i)*rpolarityinv(m,j,i) + vecp(j,i) = vecp(j,i) + + & conjp(m,i)*rpolarityinv(m,j,i) + vecs(j,i) = vecs(j,i) + + & conjs(m,i)*rpolarityinv(m,j,i) + vecps(j,i) = vecps(j,i) + + & conjps(m,i)*rpolarityinv(m,j,i) + enddo end do end if end do @@ -3746,10 +3808,20 @@ subroutine induce0d do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - zrsd(j,i) = rsd(j,i) * poli(i) - zrsdp(j,i) = rsdp(j,i) * poli(i) - zrsds(j,i) = rsds(j,i) * poli(i) - zrsdps(j,i) = rsdps(j,i) * poli(i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + zrsds(j,i) = 0.0d0 + zrsdps(j,i) = 0.0d0 + do m = 1, 3 + zrsd(j,i) = zrsd(j,i) + + & rsd(m,i)*rpolarity(m,j,i) + zrsdp(j,i) = zrsdp(j,i) + + & rsdp(m,i)*rpolarity(m,j,i) + zrsds(j,i) = zrsds(j,i) + + & rsds(m,i)*rpolarity(m,j,i) + zrsdps(j,i) = zrsdps(j,i) + + & rsdps(m,i)*rpolarity(m,j,i) + end do b = b + rsd(j,i)*zrsd(j,i) bp = bp + rsdp(j,i)*zrsdp(j,i) bs = bs + rsds(j,i)*zrsds(j,i) @@ -3805,10 +3877,16 @@ subroutine induce0d do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uind(j,i) = uind(j,i) + poli(i)*rsd(j,i) - uinp(j,i) = uinp(j,i) + poli(i)*rsdp(j,i) - uinds(j,i) = uinds(j,i) + poli(i)*rsds(j,i) - uinps(j,i) = uinps(j,i) + poli(i)*rsdps(j,i) + do m = 1, 3 + uind(j,i) = uind(j,i) + + & rpolarity(m,j,i)*rsd(m,i) + uinp(j,i) = uinp(j,i) + + & rpolarity(m,j,i)*rsdp(m,i) + uinds(j,i) = uinds(j,i) + + & rpolarity(m,j,i)*rsds(m,i) + uinps(j,i) = uinps(j,i) + + & rpolarity(m,j,i)*rsdps(m,i) + end do end do end if end do @@ -3817,7 +3895,6 @@ subroutine induce0d c c perform deallocation of some local arrays c - deallocate (poli) deallocate (rsd) deallocate (rsdp) deallocate (rsds) @@ -4708,9 +4785,8 @@ subroutine induce0e use units use uprior implicit none - integer i,j,k,iter + integer i,j,k,m,iter integer maxiter - real*8 polmin real*8 eps,epsold real*8 epsd,epsp real*8 epsds,epsps @@ -4720,7 +4796,6 @@ subroutine induce0e real*8 b,bp,bs,bps real*8 sum,sump real*8 sums,sumps - real*8, allocatable :: poli(:) real*8, allocatable :: field(:,:) real*8, allocatable :: fieldp(:,:) real*8, allocatable :: fields(:,:) @@ -4786,10 +4861,20 @@ subroutine induce0e do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - udir(j,i) = polarity(i) * field(j,i) - udirp(j,i) = polarity(i) * fieldp(j,i) - udirs(j,i) = polarity(i) * fields(j,i) - udirps(j,i) = polarity(i) * fieldps(j,i) + udir(j,i) = 0.0d0 + udirp(j,i) = 0.0d0 + udirs(j,i) = 0.0d0 + udirps(j,i) = 0.0d0 + do m = 1, 3 + udir(j,i) = udir(j,i) + + & rpolarity(m,j,i) * field(m,i) + udirp(j,i) = udirp(j,i) + + & rpolarity(m,j,i) * fieldp(m,i) + udirs(j,i) = udirs(j,i) + + & rpolarity(m,j,i) * fields(m,i) + udirps(j,i) = udirps(j,i) + + & rpolarity(m,j,i) * fieldps(m,i) + end do uind(j,i) = udir(j,i) uinp(j,i) = udirp(j,i) uinds(j,i) = udirs(j,i) @@ -4816,10 +4901,20 @@ subroutine induce0e do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uopt(k,j,i) = polarity(i) * field(j,i) - uoptp(k,j,i) = polarity(i) * fieldp(j,i) - uopts(k,j,i) = polarity(i) * fields(j,i) - uoptps(k,j,i) = polarity(i) * fieldps(j,i) + uopt(k,j,i) = 0.0d0 + uoptp(k,j,i) = 0.0d0 + uopts(k,j,i) = 0.0d0 + uoptps(k,j,i) = 0.0d0 + do m = 1, 3 + uopt(k,j,i) = uopt(k,j,i) + + & rpolarity(m,j,i) * field(m,i) + uoptp(k,j,i) = uoptp(k,j,i) + + & rpolarity(m,j,i) * fieldp(m,i) + uopts(k,j,i) = uopts(k,j,i) + + & rpolarity(m,j,i) * fields(m,i) + uoptps(k,j,i) = uoptps(k,j,i) + + & rpolarity(m,j,i) * fieldps(m,i) + end do uind(j,i) = uopt(k,j,i) uinp(j,i) = uoptp(k,j,i) uinds(j,i) = uopts(k,j,i) @@ -4868,7 +4963,6 @@ subroutine induce0e done = .false. maxiter = 100 iter = 0 - polmin = 0.00000001d0 eps = 100.0d0 c c estimated induced dipoles from polynomial predictor @@ -4896,7 +4990,6 @@ subroutine induce0e c c perform dynamic allocation of some local arrays c - allocate (poli(npole)) allocate (rsd(3,npole)) allocate (rsdp(3,npole)) allocate (rsds(3,npole)) @@ -4919,20 +5012,33 @@ subroutine induce0e call ufield0e (field,fieldp,fields,fieldps) do i = 1, npole if (douind(ipole(i))) then - poli(i) = max(polmin,polarity(i)) do j = 1, 3 - rsd(j,i) = (udir(j,i)-uind(j,i))/poli(i) - & + field(j,i) - rsdp(j,i) = (udirp(j,i)-uinp(j,i))/poli(i) - & + fieldp(j,i) - rsds(j,i) = (udirs(j,i)-uinds(j,i))/poli(i) - & + fields(j,i) - rsdps(j,i) = (udirps(j,i)-uinps(j,i))/poli(i) - & + fieldps(j,i) - zrsd(j,i) = rsd(j,i) * poli(i) - zrsdp(j,i) = rsdp(j,i) * poli(i) - zrsds(j,i) = rsds(j,i) * poli(i) - zrsdps(j,i) = rsdps(j,i) * poli(i) + rsd(j,i) = field(j,i) + rsdp(j,i) = fieldp(j,i) + rsds(j,i) = fields(j,i) + rsdps(j,i) = fieldps(j,i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + zrsds(j,i) = 0.0d0 + zrsdps(j,i) = 0.0d0 + do m = 1, 3 + rsd(j,i) = rsd(j,i) + + & (udir(m,i)-uind(m,i))*rpolarityinv(m,j,i) + rsdp(j,i) = rsdp(j,i) + + & (udirp(m,i)-uinp(m,i))*rpolarityinv(m,j,i) + rsds(j,i) = rsds(j,i) + + & (udirs(m,i)-uinds(m,i))*rpolarityinv(m,j,i) + rsdps(j,i) = rsdps(j,i) + + & (udirps(m,i)-uinps(m,i))*rpolarityinv(m,j,i) + zrsd(j,i) = zrsd(j,i) + + & rsd(m,i) * rpolarity(m,j,i) + zrsdp(j,i) = zrsdp(j,i) + + & rsdp(m,i) * rpolarity(m,j,i) + zrsds(j,i) = zrsds(j,i) + + & rsds(m,i) * rpolarity(m,j,i) + zrsdps(j,i) = zrsdps(j,i) + + & rsdps(m,i) * rpolarity(m,j,i) + end do conj(j,i) = zrsd(j,i) conjp(j,i) = zrsdp(j,i) conjs(j,i) = zrsds(j,i) @@ -4967,10 +5073,20 @@ subroutine induce0e uinp(j,i) = vecp(j,i) uinds(j,i) = vecs(j,i) uinps(j,i) = vecps(j,i) - vec(j,i) = conj(j,i)/poli(i) - field(j,i) - vecp(j,i) = conjp(j,i)/poli(i) - fieldp(j,i) - vecs(j,i) = conjs(j,i)/poli(i) - fields(j,i) - vecps(j,i) = conjps(j,i)/poli(i) - fieldps(j,i) + vec(j,i) = -field(j,i) + vecp(j,i) = -fieldp(j,i) + vecs(j,i) = -fields(j,i) + vecps(j,i) = -fieldps(j,i) + do m = 1, 3 + vec(j,i) = vec(j,i) + + & conj(m,i)*rpolarityinv(m,j,i) + vecp(j,i) = vecp(j,i) + + & conjp(m,i)*rpolarityinv(m,j,i) + vecs(j,i) = vecs(j,i) + + & conjs(m,i)*rpolarityinv(m,j,i) + vecps(j,i) = vecps(j,i) + + & conjps(m,i)*rpolarityinv(m,j,i) + end do end do end if end do @@ -5021,10 +5137,20 @@ subroutine induce0e do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - zrsd(j,i) = rsd(j,i) * poli(i) - zrsdp(j,i) = rsdp(j,i) * poli(i) - zrsds(j,i) = rsds(j,i) * poli(i) - zrsdps(j,i) = rsdps(j,i) * poli(i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + zrsds(j,i) = 0.0d0 + zrsdps(j,i) = 0.0d0 + do m = 1, 3 + zrsd(j,i) = zrsd(j,i) + + & rsd(m,i) * rpolarity(m,j,i) + zrsdp(j,i) = zrsdp(j,i) + + & rsdp(m,i) * rpolarity(m,j,i) + zrsds(j,i) = zrsds(j,i) + + & rsds(m,i) * rpolarity(m,j,i) + zrsdps(j,i) = zrsdps(j,i) + + & rsdps(m,i) * rpolarity(m,j,i) + end do b = b + rsd(j,i)*zrsd(j,i) bp = bp + rsdp(j,i)*zrsdp(j,i) bs = bs + rsds(j,i)*zrsds(j,i) @@ -5080,10 +5206,16 @@ subroutine induce0e do i = 1, npole if (douind(ipole(i))) then do j = 1, 3 - uind(j,i) = uind(j,i) + poli(i)*rsd(j,i) - uinp(j,i) = uinp(j,i) + poli(i)*rsdp(j,i) - uinds(j,i) = uinds(j,i) + poli(i)*rsds(j,i) - uinps(j,i) = uinps(j,i) + poli(i)*rsdps(j,i) + do m = 1, 3 + uind(j,i) = uind(j,i) + + & rpolarity(m,j,i)*rsd(m,i) + uinp(j,i) = uinp(j,i) + + & rpolarity(m,j,i)*rsdp(m,i) + uinds(j,i) = uinds(j,i) + + & rpolarity(m,j,i)*rsds(m,i) + uinps(j,i) = uinps(j,i) + + & rpolarity(m,j,i)*rsdps(m,i) + end do end do end if end do @@ -5092,7 +5224,6 @@ subroutine induce0e c c perform deallocation of some local arrays c - deallocate (poli) deallocate (rsd) deallocate (rsdp) deallocate (rsds) @@ -5796,7 +5927,6 @@ subroutine uscale0a (mode,rsd,rsdp,zrsd,zrsdp) real*8 xr,yr,zr real*8 r,r2,rr3,rr5 real*8 pdi,pti - real*8 polmin real*8 poli,polik real*8 damp,expdamp real*8 pgamma,off2 @@ -5817,12 +5947,16 @@ subroutine uscale0a (mode,rsd,rsdp,zrsd,zrsdp) c c use diagonal preconditioner elements as first approximation c - polmin = 0.00000001d0 do i = 1, npole - poli = udiag * max(polmin,polarity(i)) do j = 1, 3 - zrsd(j,i) = poli * rsd(j,i) - zrsdp(j,i) = poli * rsdp(j,i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + do m = 1, 3 + zrsd(j,i) = zrsd(j,i) + + & rpolarity(m,j,i)*rsd(m,i) * udiag + zrsdp(j,i) = zrsdp(j,i) + + & rpolarity(m,j,i)*rsdp(m,i) * udiag + enddo end do end do c @@ -5900,7 +6034,7 @@ subroutine uscale0a (mode,rsd,rsdp,zrsd,zrsdp) zi = z(ii) pdi = pdamp(i) pti = thole(i) - poli = polarity(i) + poli = (polarity(1,i)+polarity(2,i)+polarity(3,i))/3.0d0 do j = i+1, npole uscale(ipole(j)) = 1.0d0 end do @@ -5937,7 +6071,9 @@ subroutine uscale0a (mode,rsd,rsdp,zrsd,zrsdp) scale5 = scale5 * (1.0d0-expdamp*(1.0d0-damp)) end if end if - polik = poli * polarity(k) + polik = poli*(polarity(1,k)+ + & polarity(2,k)+ + & polarity(3,k))/3.0d0 rr3 = scale3 * polik / (r*r2) rr5 = 3.0d0 * scale5 * polik / (r*r2*r2) minv(m+1) = rr5*xr*xr - rr3 @@ -6001,7 +6137,6 @@ subroutine uscale0b (mode,rsd,rsdp,zrsd,zrsdp) real*8 xr,yr,zr real*8 r,r2,rr3,rr5 real*8 pdi,pti - real*8 polmin real*8 poli,polik real*8 damp,expdamp real*8 pgamma @@ -6029,12 +6164,16 @@ subroutine uscale0b (mode,rsd,rsdp,zrsd,zrsdp) c c use diagonal preconditioner elements as first approximation c - polmin = 0.00000001d0 do i = 1, npole - poli = udiag * max(polmin,polarity(i)) do j = 1, 3 - zrsd(j,i) = poli * rsd(j,i) - zrsdp(j,i) = poli * rsdp(j,i) + zrsd(j,i) = 0.0d0 + zrsdp(j,i) = 0.0d0 + do m = 1, 3 + zrsd(j,i) = zrsd(j,i) + + & rpolarity(m,j,i)*rsd(m,i) * udiag + zrsdp(j,i) = zrsdp(j,i) + + & rpolarity(m,j,i)*rsdp(m,i) * udiag + enddo zrsdt(j,i) = 0.0d0 zrsdtp(j,i) = 0.0d0 end do @@ -6137,7 +6276,7 @@ subroutine uscale0b (mode,rsd,rsdp,zrsd,zrsdp) zi = z(ii) pdi = pdamp(i) pti = thole(i) - poli = polarity(i) + poli = (polarity(1,i)+polarity(2,i)+polarity(3,i))/3.0d0 do j = 1, np11(ii) uscale(ip11(j,ii)) = u1scale end do @@ -6172,7 +6311,9 @@ subroutine uscale0b (mode,rsd,rsdp,zrsd,zrsdp) scale5 = scale5 * (1.0d0-expdamp*(1.0d0-damp)) end if end if - polik = poli * polarity(k) + !ACS NOTE TO JWP: the polarizability is taken as the mean here + polik = poli*(polarity(1,k)+polarity(2,k)+polarity(3,k)) + & / 3.0d0 rr3 = scale3 * polik / (r*r2) rr5 = 3.0d0 * scale5 * polik / (r*r2*r2) minv(m+1) = rr5*xr*xr - rr3 diff --git a/source/initial.f b/source/initial.f index ff6ed5490..49983cf0a 100644 --- a/source/initial.f +++ b/source/initial.f @@ -14,12 +14,6 @@ c c "initial" sets up original values for some parameters and c variables that might not otherwise get initialized -c -c note calls below to the "kmp_set" routines are for use with -c the Intel compiler, but must be commented for other compilers; -c alternatively, these values can be set via the KMP_STACKSIZE -c and KMP_BLOCKTIME environment variables -c c subroutine initial use sizes @@ -48,7 +42,6 @@ subroutine initial use scales use sequen use socket - use virial use warp use zclose implicit none @@ -82,11 +75,13 @@ subroutine initial !$ call omp_set_num_threads (nthread) !$ call omp_set_nested (.true.) c -c Intel compiler extensions to OpenMP standard, 268435456 bytes is -c 2**28 bytes, or 256 MB; comment these lines for other compilers +c Intel compiler extensions to OpenMP standard c -c!$ call kmp_set_stacksize_s (268435456) -c!$ call kmp_set_blocktime (0) +c commented out for gfortran compiler +c#ifdef __INTEL_COMPILER +c!$ call kmp_set_stacksize_s (2**28) +c!$ call kmp_set_blocktime (0) +c#endif c c values of machine precision constants c @@ -151,16 +146,12 @@ subroutine initial c use_group = .false. c -c flags for use of periodic boundaries +c flags for periodic boundaries c use_bounds = .false. use_replica = .false. use_polymer = .false. c -c flag for use of internal virial -c - use_virial = .true. -c c default values for unitcell dimensions c xbox = 0.0d0 diff --git a/source/initprm.f b/source/initprm.f index ed1fd68a0..a269edf41 100644 --- a/source/initprm.f +++ b/source/initprm.f @@ -189,7 +189,7 @@ subroutine initprm if (.not. allocated(eps4)) allocate (eps4(maxtyp)) if (.not. allocated(reduct)) allocate (reduct(maxtyp)) if (.not. allocated(chg)) allocate (chg(maxtyp)) - if (.not. allocated(polr)) allocate (polr(maxtyp)) + if (.not. allocated(polr)) allocate (polr(3,maxtyp)) if (.not. allocated(athl)) allocate (athl(maxtyp)) if (.not. allocated(pgrp)) allocate (pgrp(maxval,maxtyp)) c @@ -209,7 +209,9 @@ subroutine initprm eps4(i) = 0.0d0 reduct(i) = 0.0d0 chg(i) = 0.0d0 - polr(i) = 0.0d0 + do j = 1, 3 + polr(j,i) = 0.0d0 + enddo athl(i) = 0.0d0 do j = 1, maxval pgrp(j,i) = 0 diff --git a/source/kpolar.f b/source/kpolar.f index 18d7f7321..5a2935bef 100644 --- a/source/kpolar.f +++ b/source/kpolar.f @@ -37,13 +37,16 @@ subroutine kpolar use potent use usolve implicit none - integer i,j,k,next + integer i,j,k,kx,ky,kz,next integer nlist,npg integer pg(maxval) integer, allocatable :: list(:) - real*8 pol,thl + real*8 pol,apol(3),thl real*8 sixth logical header + character*8 axt + character*4 pa,pb,pc,pd + character*16 pt character*20 keyword character*240 record character*240 string @@ -133,6 +136,9 @@ subroutine kpolar c perform dynamic allocation of some global arrays c if (allocated(polarity)) deallocate (polarity) + if (allocated(rpolarity)) deallocate (rpolarity) + if (allocated(rpolarityinv)) deallocate (rpolarityinv) + if (allocated(is_aniso)) deallocate (is_aniso) if (allocated(thole)) deallocate (thole) if (allocated(pdamp)) deallocate (pdamp) if (allocated(udir)) deallocate (udir) @@ -140,7 +146,10 @@ subroutine kpolar if (allocated(uind)) deallocate (uind) if (allocated(uinp)) deallocate (uinp) if (allocated(douind)) deallocate (douind) - allocate (polarity(n)) + allocate (polarity(3,n)) + allocate (rpolarity(3,3,n)) + allocate (rpolarityinv(3,3,n)) + allocate (is_aniso(n)) allocate (thole(n)) allocate (pdamp(n)) allocate (udir(3,n)) @@ -199,6 +208,9 @@ subroutine kpolar record = keyline(i) call gettext (record,keyword,next) call upcase (keyword) +c +c Look for isotropic polarizabilities +c if (keyword(1:9) .eq. 'POLARIZE ') then k = 0 pol = 0.0d0 @@ -220,7 +232,9 @@ subroutine kpolar & 'Damp',5x,'Group Atom Types'/) end if if (k .le. maxtyp) then - polr(k) = pol + do j = 1, 3 + polr(j,k) = pol + enddo athl(k) = thl do j = 1, maxval pgrp(j,k) = pg(j) @@ -242,12 +256,70 @@ subroutine kpolar end if end if end if +c +c Look for anisotropic polarizabilities +c + if (keyword(1:9) .eq. 'APOLARIZE') then + k = 0 + do j = 1, 3 + apol(j) = 0.0d0 + enddo + thl = -1.0d0 + do j = 1, maxval + pg(j) = 0 + end do + call getnumb (record,k,next) + string = record(next:240) + read (string,*,err=80,end=80) apol(1),apol(2),apol(3), + & thl,(pg(j),j=1,maxval) + 80 continue + if (k .gt. 0) then + if (header .and. .not.silent) then + header = .false. + write (iout,90) + 90 format (/,' Additional Anisotropic Atomic Dipole', + & ' Polarizability Parameters :', + & //,5x,'Atom Type',9x,'Alpha xx',2x, + & 'Alpha yy',2x, 'Alpha zz',6x, + & 'Damp',4x,'Group Atom Types'/) + end if + if (k .le. maxtyp) then + do j = 1, 3 + polr(j,k) = apol(j) + enddo + athl(k) = thl + do j = 1, maxval + pgrp(j,k) = pg(j) + if (pg(j) .eq. 0) then + npg = j - 1 + goto 100 + end if + end do + 100 continue + if (.not. silent) then + write (iout,110) k,apol(1:3),thl,(pg(j),j=1,npg) + 110 format (4x,i6,10x,3f10.3,2x,f10.3,7x,20i5) + end if + else + write (iout,79) + 79 format (/,' KPOLAR -- Too many Dipole', + & ' Polarizability Parameters') + abort = .true. + end if + end if + end if end do c c find and store the atomic dipole polarizability parameters c do i = 1, n - polarity(i) = polr(type(i)) + is_aniso(i) = .false. + do j = 1, 3 + polarity(j,i) = polr(j, type(i)) + enddo + if((polarity(1,i) .ne. polarity(2,i)) .or. + & (polarity(1,i) .ne. polarity(3,i)) .or. + & (polarity(2,i) .ne. polarity(3,i))) is_aniso(i) = .true. thole(i) = athl(type(i)) end do c @@ -259,6 +331,9 @@ subroutine kpolar record = keyline(i) call gettext (record,keyword,next) call upcase (keyword) +c +c Look for isotropic polarizabilities +c if (keyword(1:9) .eq. 'POLARIZE ') then k = 0 pol = 0.0d0 @@ -267,24 +342,61 @@ subroutine kpolar if (k.lt.0 .and. k.ge.-n) then k = -k string = record(next:240) - read (string,*,err=80,end=80) pol,thl - 80 continue + read (string,*,err=120,end=120) pol,thl + 120 continue if (header) then header = .false. - write (iout,90) - 90 format (/,' Additional Dipole Polarizabilities', + write (iout,130) + 130 format (/,' Additional Dipole Polarizabilities', & ' for Specific Atoms :', & //,6x,'Atom',15x,'Alpha',8x,'Damp',/) end if if (.not. silent) then - write (iout,100) k,pol,thl - 100 format (4x,i6,10x,f10.3,2x,f10.3) + write (iout,140) k,pol,thl + 140 format (4x,i6,10x,f10.3,2x,f10.3) end if - polarity(k) = pol + do j = 1, 3 + polarity(j,k) = pol + end do + thole(k) = thl + end if + end if +c +c Look for anisotropic polarizabilities +c + if (keyword(1:9) .eq. 'APOLARIZE') then + k = 0 + do j = 1, 3 + apol(j) = 0.0d0 + enddo + thl = -1.0d0 + call getnumb (record,k,next) + if (k.lt.0 .and. k.ge.-n) then + k = -k + string = record(next:240) + read (string,*,err=150,end=150) apol(1),apol(2),apol(3), + & thl + 150 continue + if (header .and. .not.silent) then + header = .false. + write (iout,160) + 160 format (/,' Additional Anisotropic Atomic Dipole', + & ' Polarizability Parameters for Specific Atoms:', + & //,5x,'Atom Type',9x,'Alpha xx',2x, + & 'Alpha yy',2x,'Alpha zz',6x,'Damp'/) + end if + if (.not. silent) then + write (iout,170) k,apol(1:3),thl + 170 format (4x,i6,10x,3f10.3,2x,f10.3) + end if + do j = 1, 3 + polarity(j,k) = apol(j) + end do thole(k) = thl end if end if end do + c c remove zero and undefined polarizable sites from the list c @@ -292,7 +404,10 @@ subroutine kpolar if (use_polar) then npole = 0 do i = 1, n - if (polsiz(i).ne.0 .or. polarity(i).ne.0.0d0) then + if (polsiz(i) .ne. 0 .or. + & polarity(1,i) .ne. 0.0d0 .or. + & polarity(2,i) .ne. 0.0d0 .or. + & polarity(3,i) .ne. 0.0d0) then npole = npole + 1 ipole(npole) = i pollist(i) = npole @@ -303,8 +418,12 @@ subroutine kpolar do k = 1, maxpole pole(k,npole) = pole(k,i) end do - if (polarity(i) .ne. 0.0d0) npolar = npolar + 1 - polarity(npole) = polarity(i) + if (polarity(1,i).ne.0.0d0 .or. + & polarity(2,i).ne.0.0d0 .or. + & polarity(3,i).ne.0.0d0) npolar = npolar + 1 + do j = 1, 3 + polarity(j,npole) = polarity(j,i) + enddo thole(npole) = thole(i) end if end do @@ -317,7 +436,8 @@ subroutine kpolar if (thole(i) .eq. 0.0d0) then pdamp(i) = 0.0d0 else - pdamp(i) = polarity(i)**sixth + pdamp(i) = ((polarity(1,i)+polarity(2,i)+ + & polarity(3,i))/3.0d0)**sixth end if end do c diff --git a/source/kpolr.f b/source/kpolr.f index 5328a7ea6..47f538cee 100644 --- a/source/kpolr.f +++ b/source/kpolr.f @@ -20,7 +20,7 @@ module kpolr implicit none integer, allocatable :: pgrp(:,:) - real*8, allocatable :: polr(:) + real*8, allocatable :: polr(:,:) real*8, allocatable :: athl(:) save end diff --git a/source/mdsave.f b/source/mdsave.f index 3a97fa53c..e3c68ba3d 100644 --- a/source/mdsave.f +++ b/source/mdsave.f @@ -242,7 +242,9 @@ subroutine mdsave (istep,dt,epot,eksum) write (iind,280) n,title(1:ltitle) 280 format (i6,2x,a) do i = 1, npole - if (polarity(i) .ne. 0.0d0) then + if (polarity(1,i) .ne. 0.0d0 .or. + & polarity(2,i) .ne. 0.0d0 .or. + & polarity(3,i) .ne. 0.0d0) then k = ipole(i) write (iind,290) k,name(k),(debye*uind(j,i),j=1,3) 290 format (i6,2x,a3,3f12.6) diff --git a/source/mutate.f b/source/mutate.f index 5d5e5d10f..7c85f2b2e 100644 --- a/source/mutate.f +++ b/source/mutate.f @@ -190,7 +190,9 @@ subroutine altelec end do do i = 1, npolar if (mut(i)) then - polarity(i) = polarity(i) * elambda + do j = 1, 3 + polarity(j,i) = polarity(j,i) * elambda + end do end if end do end if diff --git a/source/polar.f b/source/polar.f index ce828274d..1cf532044 100644 --- a/source/polar.f +++ b/source/polar.f @@ -19,7 +19,10 @@ c optlevel current OPT order for reciprocal potential and field c copt coefficients for OPT total induced dipole moments c copm coefficients for OPT incremental induced dipole moments -c polarity dipole polarizability for each multipole site (Ang**3) +c polarity anisotropic dipole polarizability for each multipole site (body frame, Ang**3) +c rpolarity anisotropic dipole polarizability for each multipole site (lab frame, Ang**3) +c rpolarityinv inverted anisotropic dipole polarizability for each multipole site (lab frame, Ang**3) +c is_aniso whether each polarizabile site is anisotropic or not c thole Thole polarizability damping value for each site c pdamp value of polarizability scale factor for each site c udir direct induced dipole components at each multipole site @@ -49,7 +52,10 @@ module polar integer optlevel real*8, allocatable :: copt(:) real*8, allocatable :: copm(:) - real*8, allocatable :: polarity(:) + real*8, allocatable :: rpolarity(:,:,:) + real*8, allocatable :: rpolarityinv(:,:,:) + real*8, allocatable :: polarity(:,:) + logical, allocatable :: is_aniso(:) real*8, allocatable :: thole(:) real*8, allocatable :: pdamp(:) real*8, allocatable :: udir(:,:) diff --git a/source/polarize.f b/source/polarize.f index 255314a38..2579a634b 100644 --- a/source/polarize.f +++ b/source/polarize.f @@ -60,7 +60,7 @@ program polarize end if addu = 0.0d0 do i = 1, npole - addu = polarity(i) + addu + addu = (polarity(1,i)+polarity(2,i)+polarity(3,i))/3.0d0 + addu end do fstr = ' Additive Total Polarizability : ' if (nmol .eq. 1) fstr = ' Additive Molecular Polarizability :' @@ -190,7 +190,7 @@ subroutine moluind (exfield,umol) real*8 a,b,sum real*8 umol(3) real*8 exfield(3) - real*8, allocatable :: poli(:) + real*8, allocatable :: poli(:,:) real*8, allocatable :: field(:,:) real*8, allocatable :: rsd(:,:) real*8, allocatable :: zrsd(:,:) @@ -203,13 +203,13 @@ subroutine moluind (exfield,umol) c do i = 1, npole do j = 1, 3 - uind(j,i) = polarity(i) * exfield(j) + uind(j,i) = polarity(j,i) * exfield(j) end do end do c c perform dynamic allocation of some local arrays c - allocate (poli(npole)) + allocate (poli(3,npole)) allocate (field(3,npole)) allocate (rsd(3,npole)) allocate (zrsd(3,npole)) @@ -226,10 +226,10 @@ subroutine moluind (exfield,umol) polmin = 0.00000001d0 call ufield (field) do i = 1, npole - poli(i) = max(polmin,polarity(i)) do j = 1, 3 + poli(j,i) = max(polmin,polarity(j,i)) rsd(j,i) = field(j,i) - zrsd(j,i) = rsd(j,i) * poli(i) + zrsd(j,i) = rsd(j,i) * poli(j,i) conj(j,i) = zrsd(j,i) end do end do @@ -248,7 +248,7 @@ subroutine moluind (exfield,umol) do i = 1, npole do j = 1, 3 uind(j,i) = vec(j,i) - vec(j,i) = conj(j,i)/poli(i) - field(j,i) + vec(j,i) = conj(j,i)/poli(j,i) - field(j,i) end do end do a = 0.0d0 @@ -269,7 +269,7 @@ subroutine moluind (exfield,umol) b = 0.0d0 do i = 1, npole do j = 1, 3 - zrsd(j,i) = rsd(j,i) * poli(i) + zrsd(j,i) = rsd(j,i) * poli(j,i) b = b + rsd(j,i)*zrsd(j,i) end do end do diff --git a/source/poledit.f b/source/poledit.f index 9cd118839..91a56bccf 100644 --- a/source/poledit.f +++ b/source/poledit.f @@ -506,7 +506,7 @@ subroutine molsetup use mpole use polar implicit none - integer i,j,k,m,ixyz + integer i,j,jj,k,m,ixyz integer atmnum,size integer freeunit real*8 xi,yi,zi @@ -593,7 +593,7 @@ subroutine molsetup c c perform dynamic allocation of some global arrays c - if (.not. allocated(polarity)) allocate (polarity(n)) + if (.not. allocated(polarity)) allocate (polarity(3,n)) if (.not. allocated(thole)) allocate (thole(n)) if (.not. allocated(pdamp)) allocate (pdamp(n)) if (.not. allocated(ipole)) allocate (ipole(n)) @@ -605,43 +605,67 @@ subroutine molsetup do i = 1, n atmnum = atomic(i) mass(i) = 1.0d0 - polarity(i) = 0.0d0 + do j = 1, 3 + polarity(j,i) = 0.0d0 + end do thole(i) = 0.39d0 if (atmnum .eq. 1) then mass(i) = 1.008d0 - polarity(i) = 0.496d0 + do j = 1, 3 + polarity(j,i) = 0.496d0 + end do else if (atmnum .eq. 5) then mass(i) = 10.810d0 - polarity(i) = 1.600d0 + do j = 1, 3 + polarity(j,i) = 1.600d0 + end do else if (atmnum .eq. 6) then mass(i) = 12.011d0 - polarity(i) = 1.334d0 + do j = 1, 3 + polarity(j,i) = 1.334d0 + end do else if (atmnum .eq. 7) then mass(i) = 14.007d0 - polarity(i) = 1.073d0 + do j = 1, 3 + polarity(j,i) = 1.073d0 + end do else if (atmnum .eq. 8) then mass(i) = 15.999d0 - polarity(i) = 0.837d0 + do j = 1, 3 + polarity(j,i) = 0.837d0 + end do else if (atmnum .eq. 9) then mass(i) = 18.998d0 - polarity(i) = 0.507d0 + do j = 1, 3 + polarity(j,i) = 0.507d0 + end do else if (atmnum .eq. 14) then mass(i) = 28.086d0 else if (atmnum .eq. 15) then mass(i) = 30.974d0 - polarity(i) = 1.828d0 + do j = 1, 3 + polarity(j,i) = 1.828d0 + end do else if (atmnum .eq. 16) then mass(i) = 32.066d0 - polarity(i) = 3.300d0 + do j = 1, 3 + polarity(j,i) = 3.300d0 + end do else if (atmnum .eq. 17) then mass(i) = 35.453d0 - polarity(i) = 2.500d0 + do j = 1, 3 + polarity(j,i) = 2.500d0 + end do else if (atmnum .eq. 35) then mass(i) = 79.904d0 - polarity(i) = 3.595d0 + do j = 1, 3 + polarity(j,i) = 3.595d0 + end do else if (atmnum .eq. 53) then mass(i) = 126.904d0 - polarity(i) = 5.705d0 + do j = 1, 3 + polarity(j,i) = 5.705d0 + end do end if end do c @@ -652,21 +676,29 @@ subroutine molsetup if (atmnum .eq. 1) then j = i12(1,i) if (atomic(j).eq.6 .and. n12(j).eq.3) then - polarity(i) = 0.696d0 + do jj = 1, 3 + polarity(jj,i) = 0.696d0 + end do do k = 1, n12(j) m = i12(k,j) if (atomic(m).eq.8 .and. n12(m).eq.1) then - polarity(i) = 0.494d0 + do jj = 1, 3 + polarity(jj,i) = 0.494d0 + end do end if end do end if else if (atmnum .eq. 6) then if (n12(i) .eq. 3) then - polarity(i) = 1.75d0 + do jj = 1, 3 + polarity(jj,i) = 1.75d0 + end do do j = 1, n12(i) k = i12(j,i) if (atomic(k).eq.8 .and. n12(k).eq.1) then - polarity(i) = 1.334d0 + do jj = 1, 3 + polarity(jj,i) = 1.334d0 + end do end if end do end if @@ -685,7 +717,8 @@ subroutine molsetup if (thole(i) .eq. 0.0d0) then pdamp(i) = 0.0d0 else - pdamp(i) = polarity(i)**sixth + pdamp(i) = ((polarity(1,i)+polarity(2,i)+polarity(3,i)) + & /3.0d0)**sixth end if end do return @@ -1487,15 +1520,17 @@ subroutine setpolar write (iout,10) 10 format (/,' Atomic Polarizabilities for Multipole Sites :') write (iout,20) - 20 format (/,5x,'Atom',5x,'Name',7x,'Polarize',10x,'Thole',/) + 20 format (/,5x,'Atom',5x,'Name',7x,'Polarize xx',7x, + & 'Polarize yy',7x,'Polarize zz',7x,'Thole',/) do ii = 1, n i = pollist(ii) if (i .eq. 0) then write (iout,30) ii,name(ii) 30 format (i8,6x,a3,12x,'--',13x,'--') else - write (iout,40) ii,name(ii),polarity(i),thole(i) - 40 format (i8,6x,a3,4x,f12.4,3x,f12.4) + write (iout,40) ii,name(ii),polarity(1,i), + & polarity(2,i),polarity(3,i),thole(i) + 40 format (i8,6x,a3,4x,f12.4,3x,f12.4,3x,f12.4,3x,f12.4) end if end do c @@ -1529,7 +1564,9 @@ subroutine setpolar end if if (i .ne. 0) then change = .true. - polarity(i) = pol + do j = 1, 3 + polarity(j,i) = pol + enddo thole(i) = thl end if end do @@ -1540,15 +1577,17 @@ subroutine setpolar write (iout,90) 90 format (/,' Atomic Polarizabilities for Multipole Sites :') write (iout,100) - 100 format (/,5x,'Atom',5x,'Name',7x,'Polarize',10x,'Thole',/) + 100 format (/,5x,'Atom',5x,'Name',7x,'Polarize xx',7x, + & 'Polarize yy',7x,'Polarize zz',7x,'Thole',/) do ii = 1, n i = pollist(ii) if (i .eq. 0) then write (iout,110) ii,name(ii) 110 format (i8,6x,a3,12x,'--',13x,'--') else - write (iout,120) ii,name(ii),polarity(i),thole(i) - 120 format (i8,6x,a3,4x,f12.4,3x,f12.4) + write (iout,120) ii,name(ii),polarity(1,i), + & polarity(2,i),polarity(3,i),thole(i) + 120 format (i8,6x,a3,4x,f12.4,3x,f12.4,3x,f12.4,3x,f12.4) end if end do end if @@ -1779,7 +1818,7 @@ subroutine interpol real*8 polmin,norm real*8 a,b,sum real*8, allocatable :: pscale(:) - real*8, allocatable :: poli(:) + real*8, allocatable :: poli(:,:) real*8, allocatable :: field(:,:) real*8, allocatable :: rsd(:,:) real*8, allocatable :: zrsd(:,:) @@ -1795,7 +1834,7 @@ subroutine interpol c perform dynamic allocation of some local arrays c allocate (pscale(n)) - allocate (poli(npole)) + allocate (poli(3,npole)) allocate (field(3,npole)) allocate (rsd(3,npole)) allocate (zrsd(3,npole)) @@ -1807,7 +1846,7 @@ subroutine interpol call dfieldi (field,pscale) do i = 1, npole do j = 1, 3 - uind(j,i) = polarity(i) * field(j,i) + uind(j,i) = polarity(j,i) * field(j,i) end do end do c @@ -1829,10 +1868,10 @@ subroutine interpol polmin = 0.00000001d0 call ufieldi (field,pscale) do i = 1, npole - poli(i) = max(polmin,polarity(i)) do j = 1, 3 + poli(j,i) = max(polmin,polarity(j,i)) rsd(j,i) = field(j,i) - zrsd(j,i) = rsd(j,i) * poli(i) + zrsd(j,i) = rsd(j,i) * poli(j,i) conj(j,i) = zrsd(j,i) end do end do @@ -1851,7 +1890,7 @@ subroutine interpol do i = 1, npole do j = 1, 3 uind(j,i) = vec(j,i) - vec(j,i) = conj(j,i)/poli(i) - field(j,i) + vec(j,i) = conj(j,i)/poli(j,i) - field(j,i) end do end do a = 0.0d0 @@ -1872,7 +1911,7 @@ subroutine interpol b = 0.0d0 do i = 1, npole do j = 1, 3 - zrsd(j,i) = rsd(j,i) * poli(i) + zrsd(j,i) = rsd(j,i) * poli(j,i) b = b + rsd(j,i)*zrsd(j,i) end do end do @@ -2575,7 +2614,7 @@ subroutine prtpolar do j = 1, maxval if (pgrp(j,i) .ne. 0) k = j end do - write (ikey,180) ipole(i),polarity(i),thole(i), + write (ikey,180) ipole(i),polarity(1,i),thole(i), & (pgrp(j,i),j=1,k) 180 format ('polarize',2x,i5,5x,2f11.4,2x,20i5) end do diff --git a/source/prtprm.f b/source/prtprm.f index 44fa20aee..339e9214f 100644 --- a/source/prtprm.f +++ b/source/prtprm.f @@ -986,7 +986,9 @@ subroutine prtprm (itxt) c exist = .false. do i = 1, maxtyp - if (polr(i) .ne. 0.0d0) exist = .true. + do j = 1, 3 + if (polr(j,i) .ne. 0.0d0) exist = .true. + end do end do if (exist) then write (itxt,1350) formfeed,forcefield @@ -997,17 +999,19 @@ subroutine prtprm (itxt) & 6x,'Group Atom Types',/) k = 0 do i = 1, maxtyp - if (polr(i) .ne. 0.0d0) then + if (polr(1,i) .ne. 0.0d0 .or. + & polr(2,i) .ne. 0.0d0 .or. + & polr(3,i) .ne. 0.0d0) then k = k + 1 npg = 0 do j = 1, maxval if (pgrp(j,i) .ne. 0) npg = npg + 1 end do if (npg .eq. 0) then - write (itxt,1370) k,i,polr(i),athl(i) + write (itxt,1370) k,i,polr(1,i),athl(i) 1370 format (10x,i5,7x,i4,3x,2f10.3) else - write (itxt,1380) k,i,polr(i),athl(i), + write (itxt,1380) k,i,polr(1,i),athl(i), & (pgrp(j,i),j=1,npg) 1380 format (10x,i5,7x,i4,3x,2f10.3,4x,6i5) end if diff --git a/source/readprm.f b/source/readprm.f index 265ec6970..e172142b7 100644 --- a/source/readprm.f +++ b/source/readprm.f @@ -75,7 +75,7 @@ subroutine readprm real*8 an,pr,ds,dk real*8 vd,cg,dp,ps real*8 fc,bd,dl,el - real*8 pt,pol,thl + real*8 pt,pol,thl,apol(3) real*8 iz,rp,ss,ts real*8 abc,cba real*8 gi,alphi @@ -1175,7 +1175,35 @@ subroutine readprm & (pg(i),i=1,maxval) 450 continue if (ia .ne. 0) then - polr(ia) = pol + do j = 1, 3 + polr(j,ia) = pol + end do + athl(ia) = thl + do i = 1, maxval + pgrp(i,ia) = pg(i) + end do + end if + +c +c anisotropic atomic dipole polarizability parameters +c + else if (keyword(1:9) .eq. 'APOLARIZE') then + ia = 0 + do i = 1, 3 + apol(i) = 0.0d0 + end do + thl = 0.0d0 + do i = 1, maxval + pg(i) = 0 + end do + string = record(next:240) + read (string,*,err=455,end=455) ia,apol(1),apol(2),apol(3), + & thl,(pg(i),i=1,maxval) + 455 continue + if (ia .ne. 0) then + do j = 1, 3 + polr(j,ia) = apol(j) + end do athl(ia) = thl do i = 1, maxval pgrp(i,ia) = pg(i) diff --git a/source/rotpole.f b/source/rotpole.f index c9635b36d..2732d4086 100644 --- a/source/rotpole.f +++ b/source/rotpole.f @@ -29,6 +29,7 @@ subroutine rotpole do i = 1, npole call rotmat (i,a) call rotsite (i,a) + call rotalpha (i,a) end do return end @@ -252,6 +253,50 @@ subroutine rotmat (i,a) end c c +c ####################################################################### +c ## ## +c ## subroutine rotalpha -- rotate polarizabilities at single site ## +c ## ## +c ####################################################################### +c +c +c "rotalpha" computes the atomic polarizabilities at a specified site +c in the global coordinate frame by applying a rotation matrix +c +c + subroutine rotalpha (isite,a) + use sizes + use atoms + use mpole + use polar + implicit none + integer i,j,k + integer isite + real*8 a(3,3),val,valinv,kpol + real*8 polmin + + polmin = 0.00000001d0 +c +c rotate the polarizability to the global coordinate frame +c + + do i = 1, 3 + do j = 1, 3 + val = 0.0d0 + valinv = 0.0d0 + do k = 1, 3 + kpol = max(polmin,polarity(k,isite)) + val = val + a(i,k)*a(j,k)*kpol + valinv = valinv + a(i,k)*a(j,k)/kpol + end do + rpolarity(i,j,isite) = val + rpolarityinv(i,j,isite) = valinv + end do + end do + return + end +c +c c ################################################################ c ## ## c ## subroutine rotsite -- rotate multipoles at single site ## diff --git a/source/xtalfit.f b/source/xtalfit.f index d9ea131e3..9a8ed272d 100644 --- a/source/xtalfit.f +++ b/source/xtalfit.f @@ -549,7 +549,7 @@ subroutine xtalprm (mode,ixtal,xx) if (mode .eq. 'STORE') then do i = 1, npole if (type(ipole(i)) .eq. atom1) then - xx(j) = polarity(i) + xx(j) = polarity(1,i) goto 10 end if end do @@ -557,7 +557,7 @@ subroutine xtalprm (mode,ixtal,xx) sixth = 1.0d0 / 6.0d0 do i = 1, npole if (type(ipole(i)) .eq. atom1) then - polarity(i) = xx(j) + polarity(1,i) = xx(j) if (thole(i) .ne. 0.0d0) pdamp(i) = xx(j)**sixth end if end do