mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 06:05:44 -04:00
added erf(|Ri-Rj|r/rcut) factor to solvation energy and gradients...EJB
This commit is contained in:
Eric Bylaska
parent
996429f04b
commit
7f3ffb0fa2
1 changed files with 86 additions and 48 deletions
|
|
@ -122,8 +122,8 @@ c radii for heavy elements: 1.17*1.9
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
|
||||
|
|
@ -156,6 +156,9 @@ c
|
|||
if (.not.btdb_get(rtdb,'nwpw:born_dielec',mt_dbl,1,dielec))
|
||||
> dielec = 78.4d0
|
||||
|
||||
if (.not.btdb_get(rtdb,'nwpw:born_rcut',mt_dbl,1,rcut))
|
||||
> rcut = 0.1d0
|
||||
|
||||
nion = ion_nion()
|
||||
value = BA_alloc_get(mt_dbl,nion,'bradii',bradii(2),bradii(1))
|
||||
value = value.and.
|
||||
|
|
@ -198,6 +201,7 @@ c
|
|||
> " Chem. Phys. Lett., vol. 246, pages 122-129."
|
||||
write(luout,*)
|
||||
write(luout,'(5x,A,F11.6)') "dielectric constant = ",dielec
|
||||
write(luout,'(5x,A,F11.6)') "rcut = ",rcut
|
||||
if (born_relax) then
|
||||
write(luout,'(5x,A)') "self-consistent solvation"
|
||||
else
|
||||
|
|
@ -237,8 +241,8 @@ c
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
* **** local variables ****
|
||||
|
|
@ -283,8 +287,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
nwpw_born_on = born_on
|
||||
|
|
@ -303,8 +307,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
nwpw_born_relax = born_relax
|
||||
|
|
@ -331,8 +335,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
* **** local variables ****
|
||||
|
|
@ -376,7 +380,7 @@ c *** Clean up born stuff, if necessary ***
|
|||
Gsolv = nwpw_born_energy0(nion,
|
||||
> dbl_mb(ion_rion_ptr()),
|
||||
> dbl_mb(bradii(1)),
|
||||
> dbl_mb(qborn(1)),dielec)
|
||||
> dbl_mb(qborn(1)),dielec,rcut)
|
||||
if (oprint) then
|
||||
write(luout,*)
|
||||
write(luout,*) "Generalized Born Solvation"
|
||||
|
|
@ -391,6 +395,7 @@ c *** Clean up born stuff, if necessary ***
|
|||
write(luout,*)
|
||||
write(luout,'(2x,A,F8.2)') "Dielectric constant -eps- = ",
|
||||
> dielec
|
||||
write(luout,'(2x,A,F8.2)') "rcut = ",rcut
|
||||
write(luout,*)
|
||||
do ii=1,nion
|
||||
write(luout,101) ion_atom(ion_katm(ii)),ii,
|
||||
|
|
@ -424,8 +429,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
nwpw_born_screen = (1.0d0 - 1.0d0/dielec)
|
||||
|
|
@ -448,8 +453,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
* **** external functions ****
|
||||
|
|
@ -461,21 +466,21 @@ c *** Clean up born stuff, if necessary ***
|
|||
nwpw_born_energy = nwpw_born_energy0(ion_nion_qm(),
|
||||
> dbl_mb(ion_rion_ptr()),
|
||||
> dbl_mb(bradii(1)),
|
||||
> dbl_mb(qborn(1)),dielec)
|
||||
> dbl_mb(qborn(1)),dielec,rcut)
|
||||
return
|
||||
end
|
||||
|
||||
real*8 function nwpw_born_energy0(nion,rion,bradii,q,dielec)
|
||||
real*8 function nwpw_born_energy0(nion,rion,bradii,q,dielec,rcut)
|
||||
implicit none
|
||||
integer nion
|
||||
real*8 rion(3,nion),bradii(nion),q(nion)
|
||||
real*8 dielec
|
||||
real*8 dielec,rcut
|
||||
|
||||
* **** local variables ****
|
||||
integer MASTER,taskid,np
|
||||
parameter (MASTER=0)
|
||||
integer ii,jj,itask
|
||||
real*8 Gsolv,screen,C,f,dist2
|
||||
real*8 Gsolv,screen,C,f,dist2,gg
|
||||
|
||||
call Parallel_np(np)
|
||||
call Parallel_taskid(taskid)
|
||||
|
|
@ -494,7 +499,8 @@ c *** Clean up born stuff, if necessary ***
|
|||
> + (rion(3,ii)-rion(3,jj))**2)
|
||||
C = dexp(-0.25d0*dist2/(bradii(ii)*bradii(jj)))
|
||||
f = dsqrt(dist2 + bradii(ii)*bradii(jj)*C)
|
||||
Gsolv = Gsolv - 0.5d0*screen*q(ii)*q(jj)/f
|
||||
gg = erf(dsqrt(dist2)/rcut)
|
||||
Gsolv = Gsolv - 0.5d0*screen*q(ii)*q(jj)*gg/f
|
||||
end if
|
||||
itask = mod(itask+1,np)
|
||||
end do
|
||||
|
|
@ -529,8 +535,8 @@ c end do
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
* **** external functions ****
|
||||
|
|
@ -547,7 +553,7 @@ c end do
|
|||
call nwpw_born_fion0(nion,
|
||||
> dbl_mb(ion_rion_ptr()),
|
||||
> dbl_mb(bradii(1)),
|
||||
> dbl_mb(qborn(1)),dielec,dbl_mb(ftmp(1)))
|
||||
> dbl_mb(qborn(1)),dielec,rcut,dbl_mb(ftmp(1)))
|
||||
|
||||
call DAXPY_OMP(3*nion,1.0d0,dbl_mb(ftmp(1)),1,fion,1)
|
||||
|
||||
|
|
@ -558,35 +564,66 @@ c end do
|
|||
return
|
||||
end
|
||||
|
||||
real*8 function nwpw_born_gsolv(screen,qi,qj,bi,bj,xx)
|
||||
implicit none
|
||||
real*8 screen,qi,qj,bi,bj,xx
|
||||
real*8 C,f
|
||||
* **************************************
|
||||
* * *
|
||||
* * nwpw_born_gsolv *
|
||||
* * *
|
||||
* **************************************
|
||||
*
|
||||
* Calculates gsolv = -0.5*screen*qi*qj*gg/f
|
||||
*
|
||||
* where f = sqrt(xx + bi*bj*C)
|
||||
* gg = erf(sqrt(xx)/rcut)
|
||||
* C = exp(-xx/(4*bi*bj))
|
||||
*
|
||||
|
||||
real*8 function nwpw_born_gsolv(screen,qi,qj,bi,bj,xx,rcut)
|
||||
implicit none
|
||||
real*8 screen,qi,qj,bi,bj,xx,rcut
|
||||
real*8 C,f,gg
|
||||
|
||||
gg = erf(dsqrt(xx)/rcut)
|
||||
C = dexp(-0.25d0*xx/(bi*bj))
|
||||
f = dsqrt(xx + bi*bj*C)
|
||||
nwpw_born_gsolv = -0.5d0*screen*qi*qj/f
|
||||
nwpw_born_gsolv = -0.5d0*screen*qi*qj*gg/f
|
||||
return
|
||||
end
|
||||
|
||||
real*8 function nwpw_born_dgsolv(screen,qi,qj,bi,bj,xx)
|
||||
* **************************************
|
||||
* * *
|
||||
* * nwpw_born_dgsolv *
|
||||
* * *
|
||||
* **************************************
|
||||
*
|
||||
* Calculates dgsolv/dxx
|
||||
*
|
||||
* dC/dxx = -1/(4*bi*bj)*exp(-xx/(4*bi*bj)) = -1/(4*bi*bj)*C
|
||||
* df/dxx = 0.5*1/f * (dxx/dxx + bi*bj*dC/dxx) = 0.5*1/f*(1-0.25*C)
|
||||
* dgg/dxx = 1/sqrt(xx*pi) * exp(-xx/rcut**2)/rcut
|
||||
* dsqrt(xx)/dxx = 0.5/xx
|
||||
*
|
||||
real*8 function nwpw_born_dgsolv(screen,qi,qj,bi,bj,xx,rcut)
|
||||
implicit none
|
||||
real*8 screen,qi,qj,bi,bj,xx
|
||||
real*8 C,f,gsolv
|
||||
real*8 screen,qi,qj,bi,bj,xx,rcut
|
||||
real*8 C,f,gsolv,gg,dgg
|
||||
|
||||
|
||||
gg = erf(dsqrt(xx)/rcut)
|
||||
dgg = (1.0d0/dsqrt(xx*4.0d0*datan(1.0d0)))*exp(-xx/rcut**2)/rcut
|
||||
C = dexp(-0.25d0*xx/(bi*bj))
|
||||
f = dsqrt(xx + bi*bj*C)
|
||||
gsolv = -0.5d0*screen*qi*qj/f
|
||||
gsolv = -0.5d0*screen*qi*qj*gg/f
|
||||
|
||||
nwpw_born_dgsolv = -0.5d0*gsolv*(1.0-0.25d0*C)/f**2
|
||||
> - 0.5d0*screen*qi*qj*dgg/f
|
||||
return
|
||||
end
|
||||
|
||||
subroutine nwpw_born_fion0(nion,rion,bradii,q,dielec,fion)
|
||||
subroutine nwpw_born_fion0(nion,rion,bradii,q,dielec,rcut,fion)
|
||||
implicit none
|
||||
integer nion
|
||||
real*8 rion(3,nion),bradii(nion),q(nion)
|
||||
real*8 dielec
|
||||
real*8 dielec,rcut
|
||||
real*8 fion(3,nion)
|
||||
|
||||
* **** local variables ****
|
||||
|
|
@ -616,7 +653,7 @@ c end do
|
|||
|
||||
dGsolv = nwpw_born_dgsolv(screen,q(ii),q(jj),
|
||||
> bradii(ii),bradii(jj),
|
||||
> dist2)
|
||||
> dist2,rcut)
|
||||
fion(1,ii) = fion(1,ii) - 2.0d0*dGsolv*dx
|
||||
fion(2,ii) = fion(2,ii) - 2.0d0*dGsolv*dy
|
||||
fion(3,ii) = fion(3,ii) - 2.0d0*dGsolv*dz
|
||||
|
|
@ -656,8 +693,8 @@ c end do
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
integer ion_rion_ptr
|
||||
|
|
@ -666,22 +703,22 @@ c end do
|
|||
call nwpw_born_dVdq0(nion,
|
||||
> dbl_mb(ion_rion_ptr()),
|
||||
> dbl_mb(bradii(1)),
|
||||
> q,dielec,u)
|
||||
> q,dielec,rcut,u)
|
||||
return
|
||||
end
|
||||
|
||||
subroutine nwpw_born_dVdq0(nion,rion,bradii,q,dielec,u)
|
||||
subroutine nwpw_born_dVdq0(nion,rion,bradii,q,dielec,rcut,u)
|
||||
implicit none
|
||||
integer nion
|
||||
real*8 rion(3,nion),bradii(nion),q(nion)
|
||||
real*8 dielec
|
||||
real*8 dielec,rcut
|
||||
real*8 u(nion)
|
||||
|
||||
* **** local variables ****
|
||||
integer MASTER,taskid,np
|
||||
parameter (MASTER=0)
|
||||
integer ii,jj,itask
|
||||
real*8 Gsolv,screen,C,f,dist2
|
||||
real*8 Gsolv,screen,C,f,dist2,gg
|
||||
|
||||
call Parallel_np(np)
|
||||
call Parallel_taskid(taskid)
|
||||
|
|
@ -700,11 +737,12 @@ c end do
|
|||
dist2 = ((rion(1,ii)-rion(1,jj))**2
|
||||
> + (rion(2,ii)-rion(2,jj))**2
|
||||
> + (rion(3,ii)-rion(3,jj))**2)
|
||||
gg = erf(dsqrt(dist2)/rcut)
|
||||
C = dexp(-0.25d0*dist2/(bradii(ii)*bradii(jj)))
|
||||
f = dsqrt(dist2 + bradii(ii)*bradii(jj)*C)
|
||||
u(ii) = u(ii) + 0.5d0*screen*q(jj)/f
|
||||
u(jj) = u(jj) + 0.5d0*screen*q(ii)/f
|
||||
Gsolv = Gsolv - 0.5d0*screen*q(ii)*q(jj)/f
|
||||
u(ii) = u(ii) + 0.5d0*screen*q(jj)*gg/f
|
||||
u(jj) = u(jj) + 0.5d0*screen*q(ii)*gg/f
|
||||
Gsolv = Gsolv - 0.5d0*screen*q(ii)*q(jj)*gg/f
|
||||
end if
|
||||
itask = mod(itask+1,np)
|
||||
end do
|
||||
|
|
@ -732,8 +770,8 @@ c end do
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
nwpw_born_u_ptr = uborn(1)
|
||||
|
|
@ -753,8 +791,8 @@ c end do
|
|||
logical born_on,born_relax
|
||||
integer uborn(2),qborn(2)
|
||||
integer bradii(2),vradii(2),rtdb
|
||||
real*8 dielec
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,
|
||||
real*8 dielec,rcut
|
||||
common /nwpw_born_blk/ uborn,qborn,bradii,vradii,dielec,rcut,
|
||||
> rtdb,born_on,born_relax
|
||||
|
||||
nwpw_born_q_ptr = qborn(1)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue