Skip to content

Commit

Permalink
added erf(|Ri-Rj|r/rcut) factor to solvation energy and gradients...EJB
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebylaska
ebylaska committed Mar 7, 2023
1 parent 996429f commit 7f3ffb0
Showing 1 changed file with 85 additions and 47 deletions.
132 changes: 85 additions & 47 deletions src/nwpw/nwpwlib/utilities/nwpw_born.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,8 +122,8 @@ subroutine nwpw_born_init(rtdb0)
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


Expand Down Expand Up @@ -156,6 +156,9 @@ subroutine nwpw_born_init(rtdb0)
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.
Expand Down Expand Up @@ -198,6 +201,7 @@ subroutine nwpw_born_init(rtdb0)
> " 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
Expand Down Expand Up @@ -237,8 +241,8 @@ subroutine nwpw_born_end()
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 ****
Expand Down Expand Up @@ -283,8 +287,8 @@ logical function nwpw_born_on()
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
Expand All @@ -303,8 +307,8 @@ logical function nwpw_born_relax()
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
Expand All @@ -331,8 +335,8 @@ subroutine nwpw_born_Qprint(nga,nion_qm,qgaus)
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 ****
Expand Down Expand Up @@ -376,7 +380,7 @@ subroutine nwpw_born_Qprint(nga,nion_qm,qgaus)
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"
Expand All @@ -391,6 +395,7 @@ subroutine nwpw_born_Qprint(nga,nion_qm,qgaus)
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,
Expand Down Expand Up @@ -424,8 +429,8 @@ real*8 function nwpw_born_screen()
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)
Expand All @@ -448,8 +453,8 @@ real*8 function nwpw_born_energy()
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 ****
Expand All @@ -461,21 +466,21 @@ real*8 function nwpw_born_energy()
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)
Expand All @@ -494,7 +499,8 @@ real*8 function nwpw_born_energy0(nion,rion,bradii,q,dielec)
> + (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
Expand Down Expand Up @@ -529,8 +535,8 @@ subroutine nwpw_born_fion(fion)
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 ****
Expand All @@ -547,7 +553,7 @@ subroutine nwpw_born_fion(fion)
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)

Expand All @@ -558,35 +564,66 @@ subroutine nwpw_born_fion(fion)
return
end

real*8 function nwpw_born_gsolv(screen,qi,qj,bi,bj,xx)
* **************************************
* * *
* * 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
real*8 C,f
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 ****
Expand Down Expand Up @@ -616,7 +653,7 @@ subroutine nwpw_born_fion0(nion,rion,bradii,q,dielec,fion)

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
Expand Down Expand Up @@ -656,8 +693,8 @@ subroutine nwpw_born_dVdq(nion,q,u)
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
Expand All @@ -666,22 +703,22 @@ subroutine nwpw_born_dVdq(nion,q,u)
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)
Expand All @@ -700,11 +737,12 @@ subroutine nwpw_born_dVdq0(nion,rion,bradii,q,dielec,u)
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
Expand Down Expand Up @@ -732,8 +770,8 @@ integer function nwpw_born_u_ptr()
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)
Expand All @@ -753,8 +791,8 @@ integer function nwpw_born_q_ptr()
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)
Expand Down

0 comments on commit 7f3ffb0

Please sign in to comment.