comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 14:15:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Add Lebedev grids in COSMO

Browse source
This commit is contained in:
Daniel Mejia-Rodriguez 2021-11-05 19:40:22 -07:00
parent d9f2593054
commit f3df44c58d
11 changed files with 216 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/nwdft/grid/grid_checks.F
View file

@ -62,7 +62,7 @@ c
else
iangquad=mxnangpp
endif
call xc_lebedev(dum,.false.)
call xc_lebedev(dum,.false.,wa_leb,za_leb,ntheta,nphi,iangquad)
iangquad=iag_old

3
src/nwdft/grid/grid_setang.F
View file

@ -30,7 +30,8 @@ c
if (ga_nodeid().ne.0) oprint = .false.
if (leb) then ! lebedev
call xc_lebedev(l_ang, oprint)
call xc_lebedev(l_ang, oprint, wa_leb, za_leb, ntheta, nphi,
$ iangquad)
nang_leb_pts = l_ang
else ! gauss-legendre
c

12
src/nwdft/xc/xc_lebedev.F
View file

@ -1,9 +1,9 @@
subroutine xc_lebedev(l_ang, oprint)
subroutine xc_lebedev(l_ang, oprint, wa_leb, za_leb, ntheta, nphi,
$ iangquad )
implicit none
c
C$Id$
c
#include "cdft.fh"
#include "stdio.fh"
#include "global.fh"
#include "tcgmsg.fh"
@ -18,9 +18,17 @@ c Doklady Mathematics 366, 741 (1999)
c
c integer iacc ! [input]
c
integer iac_max, mx_l_ang
parameter (iac_max = 29, mx_l_ang = 5810)
integer MAXN, i
parameter(MAXN = 1281)
logical oprint
integer iangquad
double precision wa_leb(mx_l_ang,iac_max)
double precision za_leb(3,mx_l_ang,iac_max)
double precision ntheta(*), nphi(*)
double precision gjk(2,MAXN), wleb(MAXN), fasf(6), goh(3,3,48)
integer njk(6,IAC_MAX), ivj(48,6), nvj(6),
& jkv(IAC_MAX+1), itop(3,3,6), isig(3,8)

5
src/solvation/GNUmakefile
View file

@ -6,7 +6,8 @@
mnsol.o\
solv_data.o\
cosmo_smd_output.o\
cosmo_vem_output.o hnd_cosmo_lib.o
cosmo_vem_output.o hnd_cosmo_lib.o\
cosmo_lebedev.o
OBJ_OPTIMIZE = cosmo_input.o \
cosmo_initialize.o\
@ -24,7 +25,7 @@
LIB_INCLUDES = -I../nwdft/include -I../ddscf
USES_BLAS = hnd_cosmo_lib.F
USES_BLAS = hnd_cosmo_lib.F cosmo_lebedev.F
HEADERS = ../property/prop.fh ../property/atomdata.fh cosmoP.fh cosmo_params.fh

14
src/solvation/cosmo_charges.F
View file

@ -50,6 +50,7 @@ c
integer l_efciat, k_efciat
integer l_rad, k_rad, nrad
integer l_occ, k_occ
integer l_efczz, k_efczz
integer nefc, ief, i, nat
integer g_dens(3) ! 1: up, 2: down, 3: total
integer ndens,nvirt(2)
@ -91,6 +92,8 @@ c
& call errquit('cosmo_charges malloc k_efcs failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcz',l_efcz,k_efcz))
& call errquit('cosmo_charges malloc k_efcz failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efczz',l_efczz,k_efczz))
& call errquit('cosmo_charges malloc k_efczz failed',911,MA_ERR)
if(.not.ma_push_get(mt_int,nefc,'cosmo efciat',l_efciat,k_efciat))
& call errquit('cosmo_charges malloc k_efciat failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nat,'cosmo rad',l_rad,k_rad))
@ -105,6 +108,9 @@ c
if(.not.rtdb_get(rtdb,'cosmo:efcs',mt_dbl, nefc,dbl_mb(k_efcs)))
& call errquit('cosmo_charges: rtdb get failed for efcs ',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efczz',mt_dbl,nefc,dbl_mb(k_efczz)))
& call errquit('cosmo_charges: rtdb get failed for efczz',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efciat',
& mt_int,nefc,int_mb(k_efciat)))
& call errquit('cosmo_charges: rtdb get failed for efciat',914,
@ -191,6 +197,7 @@ c ----- get -cosmo- charges -----
c
call hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& dbl_mb(k_efcc),dbl_mb(k_efcs),dbl_mb(k_efcz),
& dbl_mb(k_efczz),
& int_mb(k_efciat),dbl_mb(k_rad),ecos,cosmo_file)
c
c ----- release memory block -----
@ -250,6 +257,7 @@ c
integer l_efciat, k_efciat
integer l_rad, k_rad, nrad
integer l_occ, k_occ
integer l_efczz, k_efczz
integer nefc, ief, i, nat
integer g_dens(3) ! 1: up, 2: down, 3: total
integer ndens,nvirt(2)
@ -293,6 +301,8 @@ c
& call errquit('cosmo_charges malloc k_efcs failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcz',l_efcz,k_efcz))
& call errquit('cosmo_charges malloc k_efcz failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efczz',l_efczz,k_efczz))
& call errquit('cosmo_charges malloc k_efczz failed',911,MA_ERR)
if(.not.ma_push_get(mt_int,nefc,'cosmo efciat',l_efciat,k_efciat))
& call errquit('cosmo_charges malloc k_efciat failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nat,'cosmo rad',l_rad,k_rad))
@ -306,6 +316,9 @@ c
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efcs',mt_dbl, nefc,dbl_mb(k_efcs)))
& call errquit('cosmo_charges: rtdb get failed for efcs ',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efczz',mt_dbl, nefc,dbl_mb(k_efczz)))
& call errquit('cosmo_charges: rtdb get failed for efczz',914,
& RTDB_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efciat',
& mt_int,nefc,int_mb(k_efciat)))
@ -349,6 +362,7 @@ c ----- get -cosmo- charges -----
c
call hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& dbl_mb(k_efcc),dbl_mb(k_efcs),dbl_mb(k_efcz),
& dbl_mb(k_efczz),
& int_mb(k_efciat),dbl_mb(k_rad),ecos,cosmo_file)
c
c ----- release memory block -----

39
src/solvation/cosmo_initialize.F
View file

@ -263,6 +263,7 @@ c
gammas =1.00d+00 ! full switching, see York, Karplus Eq.86.
swtol =1.00d-04 ! switching tolerance, see York, Karplus
! Eq.62-88.
iangleb =-1
c
c MN solvation models -->
c
@ -316,6 +317,7 @@ c
status = rtdb_get(rtdb,'cosmo:gamma_s',mt_dbl,1,gammas)
status = rtdb_get(rtdb,'cosmo:sw_tol',mt_dbl,1,swtol)
status = rtdb_cget(rtdb,'cosmo:solvent',1,solvname_short)
status = rtdb_get(rtdb,'cosmo:iangleb',mt_int,1,iangleb)
c
c nrad = 0 ! this block was moved below
c if (.not. rtdb_get(rtdb,'cosmo:nrrad',mt_int,1,nrad)) nrad = 0
@ -535,6 +537,20 @@ c
endif
enddo
c
c ---- lebedev ----
c
if ((do_cosmo_model.eq.do_cosmo_ks) .and. (iangleb.gt.0) ) then
iangleb = -1
endif
if (zeta.eq.0d0) then
if (iangleb.gt.0) then
zeta = 4.907d0
else
zeta = 4.500d0
endif
endif
c
c ----- create boundary elements of cavity -----
c
call hnd_cosset(rtdb,nat,dbl_mb(k_coscoor),dbl_mb(k_rad))
@ -553,14 +569,19 @@ c
c
if(osome.and.me.eq.0) then
write(luout,9993) lineq
write(luout,9991) minbem
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(luout,9992) maxbem
endif
if(ificos.eq.0) then
write(luout,9990)
if (iangleb.lt.1) then
write(luout,9991) minbem
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(luout,9992) maxbem
endif
if(ificos.eq.0) then
write(luout,9990)
else
write(luout,9989)
endif
else
write(luout,9989)
write(luout,9969)
write(luout,9968) iangleb
endif
if (do_cosmo_model.eq.DO_COSMO_KS) then
write(luout,9994) rsolv
@ -685,6 +706,10 @@ c
9972 format(' surface charge correction = ',a)
9971 format(' Cosmo: York-Karplus, doi: 10.1021/jp992097l')
9970 format(' Cosmo: Klamt-Schuurmann, doi: 10.1039/P29930000799')
9969 format(' Using -Lebedev- grid')
9968 format(' -iangleb',I4)
end
C>
C> @}

11
src/solvation/cosmo_input.F
View file

@ -72,6 +72,7 @@ c
integer iat, nfield, nrad, irad
integer iscren, ifscrn, minbem, maxbem, ificos, lineq
integer cosmo_screen, cosmo_sccor
integer iangleb
integer l_rad, k_rad
double precision dielec,dielecinf,rsolv,rad,zeta,gammas,swtol
character*255 field
@ -140,6 +141,7 @@ c
zeta =4.50d+00 ! apparently for Lebedev grids zeta=4.9
gammas =1.00d+00 ! full switching, see York, Karplus Eq.86.
swtol =1.00d-04 ! switching tolerance, see York, Karplus Eq.62-88.
iangleb =-1
c
c MN solvation models -->
c
@ -492,6 +494,15 @@ c
if (.not. rtdb_put(rtdb,'cosmo:ificos',mt_int,1,ificos))
$ call errquit('cosmo_input: rtdb put failed',911, RTDB_ERR)
c
c --- lebedev
c
else if(inp_compare(.false.,'lebedev',field)) then
if(.not. inp_i(iangleb)) iangleb = 4
if (iangleb.lt.1)
$ call errquit('cosmo_input: lebedev <1',911,INPUT_ERR)
if(.not. rtdb_put(rtdb,'cosmo:iangleb',mt_int,1,iangleb))
$ call errquit('cosmo_input: rtdb put failed',911,RTDB_ERR)
c
c --- lineq
c
else if(inp_compare(.false.,'lineq',field)) then

35
src/solvation/cosmo_lebedev.F Normal file
View file

@ -0,0 +1,35 @@
subroutine cosmo_lebedev(iang,n,w,xyz)
implicit none
#include "errquit.fh"
integer iac_max, mx_l_ang
parameter (iac_max=29, mx_l_ang=5810)
integer iang,jang,i,n
integer idum(mx_l_ang)
double precision xyz(3,*),w(*)
double precision wleb(mx_l_ang,iac_max)
double precision zleb(3,mx_l_ang,iac_max)
integer, external :: xc_leb_acc2nrd
if (iang.lt.1)
$ call errquit('cosmo_lebedev: iangleb is <1',iang,INPUT_ERR)
if (iang.gt.iac_max)
$ call errquit('cosmo_lebedev: iangleb is >29',iang,INPUT_ERR)
select case (iang)
case (3,4)
jang = 4
case (9,10,11)
jang = 11
case default
jang = iang
end select
call xc_lebedev(n, .false., wleb, zleb, idum, idum, jang)
n = xc_leb_acc2nrd(jang)
call dcopy(3*n,zleb(1,1,jang),1,xyz,1)
call dcopy(n,wleb(1,jang),1,w,1)
end

4
src/solvation/cosmo_params.fh
View file

@ -4,10 +4,10 @@ c
double precision dielec,dielecinf,screen,rsolv,zeta,gammas,swtol
double precision adiag,dsurf,dvol,srfmol,volmol,ptspatm
integer do_cosmo_model,lineq,minbem,maxbem,ificos,ifscrn,
& cosmo_sccor
& cosmo_sccor,iangleb
common/hnd_cospar/dielec,dielecinf,screen,rsolv,zeta,gammas,swtol
common/hnd_cosmod/do_cosmo_model,lineq,minbem,maxbem,
I ificos,ifscrn,
& cosmo_sccor
& cosmo_sccor,iangleb
common/hnd_cosdat/adiag,dsurf,dvol,srfmol,volmol,ptspatm

19
src/solvation/grad_hnd_cos.F
View file

@ -110,6 +110,8 @@ C-------------------------local variables--------------------------
integer iefc ! counter over COSMO charges
integer iefc_c ! memory index for COSMO charge coordinates
integer iefc_q ! memory index for COSMO charges
integer l_efczz
integer k_efczz
integer max_sh_bf
integer l_efcs, k_efcs
@ -198,7 +200,11 @@ c
if(.not.ma_push_get(mt_dbl,nefc,'efcs',l_efcs,k_efcs))
$ call errquit('grad_hnd_cos: efcs alloc failed',0,MA_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efcs',mt_dbl,nefc,dbl_mb(k_efcs)))
$ call errquit('grad_hnd_cos: rtdb_get failed',0,RTDB_ERR)
$ call errquit('grad_hnd_cos: rtdb_get failed',0,RTDB_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'efczz',l_efczz,k_efczz))
$ call errquit('grad_hnd_cos: efczz alloca failed',0,MA_ERR)
if(.not.rtdb_get(rtdb,'cosmo:efczz',mt_dbl,nefc,dbl_mb(k_efczz)))
$ call errquit('grad_hnd_cos: rtdb_get failed',0,RTDB_ERR)
if(.not.ma_push_get(mt_dbl,max_sh_bf**2,'shell p',l_p,k_p))
$ call errquit('grad_hnd_cos: shell p alloc failed',0,MA_ERR)
if(.not.ma_push_get(mt_int,2*nat,'atom->efc map',l_map,k_map))
@ -425,17 +431,13 @@ c
invscreen = 1.0d0/(1.0d0*screen)
do ich1 = ga_nodeid()+1, nefc, nproc
if (do_cosmo_model.eq.DO_COSMO_YK) then
zeta1 = zeta*sqrt(ptspatm)
$ / (dbl_mb(k_rad+int_mb(k_efciat+ich1-1)-1)
$ * sqrt(4.0d0*pi))
zeta1 = dbl_mb(k_efczz+ich1-1)
endif
do ich2 = 1, ich1-1
if (int_mb(k_efciat+ich1-1).eq.int_mb(k_efciat+ich2-1)) cycle
if (do_cosmo_model.eq.DO_COSMO_YK) then
zeta2 = zeta*sqrt(ptspatm)
$ / (dbl_mb(k_rad+int_mb(k_efciat+ich2-1)-1)
$ * sqrt(4.0d0*pi))
zeta2 = dbl_mb(k_efczz+ich2-1)
zeta12 = zeta1*zeta2/sqrt(zeta1**2+zeta2**2)
endif
dx = dbl_mb(0+3*(ich1-1)+iefc_c)
@ -493,6 +495,9 @@ c
if (.not.ma_pop_stack(l_p))
$ call errquit("grad_hnd_cos: could not deallocate l_p",
$ 0,MA_ERR)
if (.not.ma_pop_stack(l_efczz))
$ call errquit("grad_hnd_cos: could not deallocate l_efczz",
$ 0,MA_ERR)
if (.not.ma_pop_stack(l_efcs))
$ call errquit("grad_hnd_cos: could not deallocate l_efcs",
$ 0,MA_ERR)

152
src/solvation/hnd_cosmo_lib.F
View file

@ -179,6 +179,7 @@ c
integer numpps( mxfac,mxatm)
double precision apex(3,mxapx)
double precision xyzff( mxfac,mxatm)
double precision wleb(mxfac)
c
logical some
logical out
@ -205,9 +206,15 @@ c
nseg = 0
nfac = 0
ndiv = 0
call hnd_cossph(nseg,nfac,ndiv,
1 ijkfac,xyzseg,ijkseg,mxfac,apex,mxapx,
2 dsurf,dvol,adiag)
if (iangleb.lt.1) then
call hnd_cossph(nseg,nfac,ndiv,
1 ijkfac,xyzseg,ijkseg,mxfac,apex,mxapx,
2 dsurf,dvol,adiag)
else
call cosmo_lebedev(iangleb, nseg, wleb, xyzseg)
endif
ptspatm = dble(nseg)
c
c ----- debug printing -----
@ -222,7 +229,7 @@ c
1 xyzseg(3,lseg),ijkseg( lseg)
enddo
endif
if(dbug) then
if(dbug .and. iangleb.lt.1) then
write(luout,9996)
do lfac=1,nfac
mfac=lfac+nseg
@ -257,7 +264,7 @@ c
call hnd_cossas(nat,xyzatm,ratm,mxatm,
1 nspa,nppa,xyzspa,ijkspa,
2 nseg,nfac,xyzseg,ijkseg,insseg,
3 numpps,xyzff,mxfac,rtdb)
3 numpps,xyzff,wleb,mxfac,rtdb)
c
return
9999 format(' nseg,nfac,ndiv=nfac/nseg,dsurf,dvol = ',3i7,2f10.6)
@ -350,7 +357,7 @@ C>
subroutine hnd_cossas(nat,xyzatm,ratm,mxatom,
1 nspa,nppa,xyzspa,ijkspa,
2 nseg,nfac,xyzseg,ijkseg,insseg,
3 numpps,xyzff,mxface,rtdb)
3 numpps,xyzff,wleb,mxface,rtdb)
implicit none
#include "cosmo_params.fh"
#include "errquit.fh"
@ -395,10 +402,12 @@ c
integer ijkspa( mxface,mxatom)
integer numpps( mxface,mxatom)
double precision xyzff( mxface,mxatom)
double precision wleb( mxface)
double precision zero, one
data one /1.0d+00/
integer l_efcc, k_efcc, l_efcs, k_efcs, l_efcz, k_efcz
integer l_efclb, k_efclb, k_efciat, l_efciat
integer l_efczz, k_efczz
double precision ratm_real,dij,dum,cavdsp,pi,zetai,zetaii
integer m,mfac,mseg
integer nefc,iat,jat,npp,i,iseg,ifac,ief,ipp
@ -515,7 +524,7 @@ c
endif
endif
enddo
if(dbug.and.ga_nodeid().eq.0) then
if(dbug.and.ga_nodeid().eq.0.and.iangleb.lt.1) then
write(luout,9996) iat
write(luout,9995) (ijkspa(ifac+nseg,iat),ifac=1,nfac)
endif
@ -617,6 +626,8 @@ c
& call errquit('cosmo_cossas malloc k_efciat failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcz',l_efcz,k_efcz))
& call errquit('cosmo_cossas malloc k_efcz failed',911,MA_ERR)
if(.not.ma_push_get(mt_dbl,nefc,'cosmo zeta',l_efczz,k_efczz))
& call errquit('cosmo_cossas malloc k_efczz failed',911,MA_ERR)
if(.not.ma_push_get(mt_byte,nefc*8,'cosmo tags',l_efclb,k_efclb))
& call errquit('cosmo_cossas malloc k_tag failed',911,MA_ERR)
c
@ -651,15 +662,24 @@ c
c
c --- eval eq.(67) from [2] ---
c
dum=4.00d0**(maxbem-minbem)
dum=16.0d0 ! MAXBEM is obsolete in York and Karplus approach
zetai=zeta/sqrt(dsurf)/ratm_real
zetaii=zetai*sqrt(2d0/pi)
dbl_mb(k_efcs+ief-1) = (zetaii/xyzff(iseg,iat))
srfmol = srfmol + xyzff(iseg,iat)*dsurf*ratm_real**2
volmol = volmol + xyzff(iseg,iat)*dvol *ratm_real**3
if (iangleb.lt.1) then
zetai=zeta/sqrt(dsurf)/ratm_real
zetaii=zetai*sqrt(2d0/pi)
dbl_mb(k_efcs+ief-1) = (zetaii/xyzff(iseg,iat))
srfmol = srfmol+xyzff(iseg,iat)*dsurf*ratm_real**2
volmol = volmol+xyzff(iseg,iat)*dvol *ratm_real**3
else
zetai=zeta/(ratm_real*sqrt(wleb(iseg)))
zetaii=zetai*sqrt(2d0/pi)
dbl_mb(k_efcs+ief-1) = zetaii/xyzff(iseg,iat)
srfmol = srfmol + xyzff(iseg,iat)*
$ wleb(iseg)*ratm_real**2
volmol = volmol + xyzff(iseg,iat)*
$ wleb(iseg)*ratm_real**3/3d0
endif
endif
int_mb(k_efciat+ief-1)=iat
dbl_mb(k_efczz+ief-1) =zetai
endif
enddo
endif
@ -725,6 +745,9 @@ c
if(.not.rtdb_put(rtdb,'cosmo:efcs',mt_dbl, nefc,dbl_mb(k_efcs)))
$ call errquit('hnd_cossas: rtdb put failed for efcs ',914,
& rtdb_err)
if(.not.rtdb_put(rtdb,'cosmo:efczz',mt_dbl,nefc,dbl_mb(k_efczz)))
& call errquit('hnd_cossas: rtdb put failed for efczz ',914,
& rtdb_err)
c
c ----- reset cosmo:rawt to avoid trouble in cosmo charge
c calculation -----
@ -748,6 +771,8 @@ c & rtdb_err)
c
if(.not.ma_pop_stack(l_efclb)) call
& errquit('cosmo_cossas chop stack k_efclb failed',911,MA_ERR)
if(.not.ma_pop_stack(l_efczz)) call
& errquit('cosmo_cossas chop stack k_efczz failed',911,MA_ERR)
if(.not.ma_pop_stack(l_efcz)) call
& errquit('cosmo_cossas chop stack k_efcz failed',911,MA_ERR)
if(.not.ma_pop_stack(l_efciat)) call
@ -2007,7 +2032,7 @@ C> <a href="https://doi.org/10.1111/1468-0262.00296">
C> 10.1111/1468-0262.00296</a>.
C>
subroutine hnd_coschg(g_dens,ndens,rtdb,geom,basis,nat,nefc,
& efcc,efcs,efcz,efciat,ratm,ecos,
& efcc,efcs,efcz,efczz,efciat,ratm,ecos,
& cosmo_file_in)
implicit none
#include "errquit.fh"
@ -2040,6 +2065,8 @@ c
!< surface area
c
double precision efcz(nefc) !< [Output] the COSMO charges
double precision efczz(nefc) !< [Input[ the COSMO zeta value for
!< each surface charge
integer efciat(nefc) !< [Input] the atom associated
!< with each surface charge
double precision ratm(nat) !< [Input] the atom radii
@ -2399,7 +2426,7 @@ c
call dcopy(nefc,x(iep),1,x(ieq),1)
call dfill(nefc,1d0,x(ieq+nefc),1)
call hnd_cosmata(nat,nefc,efcc,efcs,efciat,ratm,x(i30))
call hnd_cosmata(nat,nefc,efcc,efcs,efczz,efciat,ratm,x(i30))
c
nodcmp=0
call hnd_linequ(x(i30),nefc,x(ieq),nefc,
@ -2461,7 +2488,7 @@ c
call hnd_cg(nat,x(iep),x(ieq),nefc,
1 x(i40),x(i50),x(i60),x(i70),
2 efcc,efcs,efciat,ratm)
2 efcc,efcs,efczz,efciat,ratm)
c
if (all) then
if (.not.rtdb_put(rtdb,'cosmo:qraw',mt_dbl,nefc,x(i20))) then
@ -2495,7 +2522,7 @@ c
call dfill(nefc,one,x(i30),1)
call hnd_cg(nat,x(i30),x(i80),nefc,
1 x(i40),x(i50),x(i60),x(i70),
2 efcc,efcs,efciat,ratm)
2 efcc,efcs,efczz,efciat,ratm)
c
if (.not.rtdb_put(rtdb,'cosmo:rawt',mt_dbl,nefc,x(i80)))
1 call errquit('hnd_coschg: could not store raw A^{-1}1',
@ -2625,8 +2652,8 @@ c identical. If surface charge corrections are employed using
c scaling factors or Lagrangian multipliers then additional
c correction terms are needed to ensure the expressions agree.
c
call hnd_cos_energy(nat,nefc,chgfac,efcc,efcs,efcz,efciat,ratm,
& dbl_mb(k_xyzpt),dbl_mb(k_zanpt),x(i10),
call hnd_cos_energy(nat,nefc,chgfac,efcc,efcs,efcz,efczz,efciat,
& ratm,dbl_mb(k_xyzpt),dbl_mb(k_zanpt),x(i10),
& allefc,atmefc,elcefc,efcefc)
solnrg= allefc+efcefc
c COSMO contribution Alternative 1 (Ref.[1] Eq.(8))
@ -2661,7 +2688,13 @@ c
solnrg= pt5* allefc
c COSMO contribution Alternative 2 (Ref.[1] Eq.(11))
c ecos = pt5*(atmefc-elcefc)+elambda
ecos = pt5*(atmefc-elcefc)
if (do_cosmo_smd) then
ecos = pt5*(atmefc-elcefc)
else
ecos = efcefc + atmefc
endif
if (oprint_energies) then
write(luout,*)'Alternative 2'
write(luout,9991) atmefc
@ -3080,8 +3113,8 @@ C> pp 11060-11079, DOI:
C> <a href="https://doi.org/10.1021/jp992097l">
C> 10.1021/jp992097l</a>.
C>
subroutine hnd_cos_energy(nat,nefc,chgscr,efcc,efcs,efcz,efciat,
& ratm,catm,zatm,pot,allefc,atmefc,elcefc,efcefc)
subroutine hnd_cos_energy(nat,nefc,chgscr,efcc,efcs,efcz,efczz,
& efciat,ratm,catm,zatm,pot,allefc,atmefc,elcefc,efcefc)
implicit none
c
#include "cosmoP.fh"
@ -3104,6 +3137,7 @@ c
!< (see e.g. Eq.(7b)) or Ref.[2]
!< Eq.(67).
double precision efcz(nefc) !< [Input] The COSMO charges
double precision efczz(nefc)
double precision ratm(nat) !< [Input] The atomic radii
double precision catm(3,nat) !< [Input] The atomic coordinates
double precision zatm(nat) !< [Input] The nuclear charges
@ -3181,12 +3215,11 @@ c
yj=efcc(2,jef)
zj=efcc(3,jef)
qj=efcz( jef)
zetaj=zeta*sqrt(ptspatm)/(ratm(efciat(jef))*sqrt(4.0d0*pi))
zetaj=efczz(jef)
efcefc=efcefc+qj*efcs(jef)*qj
c
do ief=jef+1,nefc
zetai=zeta*sqrt(ptspatm)
& /(ratm(efciat(ief))*sqrt(4.0d0*pi))
zetai=efczz(ief)
xi=efcc(1,ief)
yi=efcc(2,ief)
zi=efcc(3,ief)
@ -3212,7 +3245,7 @@ C>
C> Compute matrix `A` as needed for the in-memory COSMO
C> charge fitting.
C>
subroutine hnd_cosmata(nat,nefc,efcc,efcs,efciat,ratm,a)
subroutine hnd_cosmata(nat,nefc,efcc,efcs,efczz,efciat,ratm,a)
implicit none
#include "cosmoP.fh"
#include "cosmo_params.fh"
@ -3220,6 +3253,7 @@ C>
integer nefc !< [Input] The number of surface charges
double precision efcc(3,nefc) !< [Input] The surface charge coords
double precision efcs(nefc) !< [Input] The surface areas
double precision efczz(nefc)
integer efciat(nat) !< [Input] The atom of a surface
!< charge
double precision ratm(nat) !< [Input] The atom radii
@ -3259,26 +3293,23 @@ c
factor = zeta*dsqrt(ptspatm/(4d0*pi))
factor2 = 2d0/sqrt(pi)
do jef=1,nefc
zetaj=factor/ratm(efciat(jef))
zetaj=efczz(jef)
xj=efcc(1,jef)
yj=efcc(2,jef)
zj=efcc(3,jef)
a(jef,jef)=efcs(jef)
do ief=1,nefc
if(ief.eq.jef) then
aii=efcs(ief)
a(ief,jef)=aii
if(ief.eq.jef) cycle
zetai=efczz(ief)
zetaij=zetai*zetaj/sqrt(zetai**2+zetaj**2)
xi=efcc(1,ief)
yi=efcc(2,ief)
zi=efcc(3,ief)
dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
if (dij.lt.1.0d-5) then
a(ief,jef)=factor2*zetaij*(1d0 - (zetaij*dij)**2/3d0)
else
zetai=factor/ratm(efciat(ief))
zetaij=zetai*zetaj/sqrt(zetai**2+zetaj**2)
xi=efcc(1,ief)
yi=efcc(2,ief)
zi=efcc(3,ief)
dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
if (dij.lt.1.0d-5) then
a(ief,jef)=factor2*zetaij*(1d0 - (zetaij*dij)**2/3d0)
else
a(ief,jef)=derf(zetaij*dij)/dij
endif
a(ief,jef)=derf(zetaij*dij)/dij
endif
enddo
enddo
@ -3293,7 +3324,7 @@ C> for the charge vector using a dot-product based algorithm. The matrix
C> is generated on the fly. For performance reasons the routine is
C> replicated data parallel.
C>
subroutine hnd_cosaxd(nat,x,ax,nefc,efcc,efcs,efciat,ratm)
subroutine hnd_cosaxd(nat,x,ax,nefc,efcc,efcs,efczz,efciat,ratm)
implicit none
#include "global.fh"
#include "msgids.fh"
@ -3304,6 +3335,7 @@ c
integer nefc !< [Input] The number of surface charges
double precision efcc(3,nefc) !< [Input] The surface charge coords
double precision efcs(nefc) !< [Input] The surface areas
double precision efczz(nefc)
integer efciat(nefc) !< [Input] The atom of a surface
!< charge
double precision ratm(nat) !< [Input] The atom radii
@ -3344,11 +3376,11 @@ c
factor = zeta*dsqrt(ptspatm/(4d0*pi))
factor2 = 2d0/sqrt(pi)
do i=ga_nodeid()+1,nefc,ga_nnodes()
zetai = ratm(efciat(i))**2
zetai = efczz(i)
ax(i) = ax(i) + efcs(i)*x(i)
do j=i+1,nefc
zetaj = ratm(efciat(j))**2
zetaij= factor/sqrt(zetai+zetaj)
zetaj = efczz(j)
zetaij= zeta*zetaj/sqrt(zetai**2+zetaj**2)
dij=r(i,j)
if (dij.lt.1.0d-5) then
aij=factor2*zetaij*(1d0 - (zetaij*dij)**2/3d0)
@ -3373,12 +3405,13 @@ C> with the matrix. The routine is replicated data parallel.
C> The matrix `A` is symmetric so we simply call the matrix-vector
C> product.
C>
subroutine hnd_cosxad(nat,x,xa,nefc,efcc,efcs,efciat,ratm)
subroutine hnd_cosxad(nat,x,xa,nefc,efcc,efcs,efczz,efciat,ratm)
implicit none
integer nat !< [Input] The number of atoms
integer nefc !< [Input] The number of surface charges
double precision efcc(3,nefc) !< [Input] The surface charge coords
double precision efcs(nefc) !< [Input] The surface areas
double precision efczz(nefc)
integer efciat(nefc) !< [Input] The atom of a surface
!< charge
double precision ratm(nat) !< [Input] The atom radii
@ -3386,7 +3419,7 @@ C>
double precision xa(nefc) !< [Output] Vector-matrix product
!< `xA`
c
call hnd_cosaxd(nat,x,xa,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,x,xa,nefc,efcc,efcs,efczz,efciat,ratm)
c
return
end
@ -3411,8 +3444,8 @@ C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI:
C> <a href="https://doi.org/10.1039/P29930000799">
C> 10.1039/P29930000799</a>.
C>
subroutine hnd_cosequ(nat,b,x,nefc,g,h,xi,xj,efcc,efcs,efciat,
& ratm)
subroutine hnd_cosequ(nat,b,x,nefc,g,h,xi,xj,efcc,efcs,efczz,
& efciat,ratm)
implicit none
c
c ----- solve A * x = b , using an iterative procedure -----
@ -3429,6 +3462,7 @@ c
integer nefc !< [Input] The number of surface charges
double precision efcc(3,nefc) !< [Input] The surface charge coords
double precision efcs(nefc) !< [Input] The surface areas
double precision efczz(nefc)
integer efciat(nefc) !< [Input] The atom of a surface
!< charge
double precision ratm(nat) !< [Input] The atom radii
@ -3458,7 +3492,7 @@ c
eps2=nefc*eps*eps
irst=0
10 irst=irst+1
call hnd_cosaxd(nat,x,xi,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,x,xi,nefc,efcc,efcs,efczz,efciat,ratm)
rp=zero
bsq=zero
do j=1,nefc
@ -3466,13 +3500,13 @@ c
xi(j)=xi(j)-b(j)
rp=rp+xi(j)*xi(j)
enddo ! j
call hnd_cosxad(nat,xi,g,nefc,efcc,efcs,efciat,ratm)
call hnd_cosxad(nat,xi,g,nefc,efcc,efcs,efczz,efciat,ratm)
do j=1,nefc
g(j)=-g(j)
h(j)= g(j)
enddo ! j
do iter=1,10*nefc
call hnd_cosaxd(nat,h,xi,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,h,xi,nefc,efcc,efcs,efczz,efciat,ratm)
anum=zero
aden=zero
do j=1,nefc
@ -3488,7 +3522,7 @@ c
xi(j)=x(j)
x(j)=x(j)+anum*h(j)
enddo ! j
call hnd_cosaxd(nat,x,xj,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,x,xj,nefc,efcc,efcs,efczz,efciat,ratm)
rsq=zero
do j=1,nefc
xj(j)=xj(j)-b(j)
@ -3509,7 +3543,7 @@ c
go to 10
endif
rp=rsq
call hnd_cosxad(nat,xj,xi,nefc,efcc,efcs,efciat,ratm)
call hnd_cosxad(nat,xj,xi,nefc,efcc,efcs,efczz,efciat,ratm)
gg=zero
dgg=zero
do j=1,nefc
@ -3660,7 +3694,8 @@ c
return
end
subroutine hnd_cg(nat,b,x,nefc,y,r,p,ap,efcc,efcs,efciat,ratm)
subroutine hnd_cg(nat,b,x,nefc,y,r,p,ap,efcc,efcs,efczz,efciat,
$ ratm)
#ifdef USE_OPENMP
use omp_lib
#endif
@ -3674,6 +3709,7 @@ c
integer nefc !< [Input] The number of surface charges
double precision efcc(3,nefc) !< [Input] The surface charge coords
double precision efcs(nefc) !< [Input] The surface areas
double precision efczz(nefc)
integer efciat(nefc) !< [Input] The atom of a surface
!< charge
double precision ratm(nat) !< [Input] The atom radii
@ -3705,7 +3741,7 @@ c
if (xnorm.eq.0d0) then
r(:) = b(:)
else
call hnd_cosaxd(nat,x,r,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,x,r,nefc,efcc,efcs,efczz,efciat,ratm)
r(:) = b(:) - r(:)
endif
@ -3724,7 +3760,7 @@ c
do while(iter.le.nefc)
iter = iter + 1
call hnd_cosaxd(nat,p,ap,nefc,efcc,efcs,efciat,ratm)
call hnd_cosaxd(nat,p,ap,nefc,efcc,efcs,efczz,efciat,ratm)
alpha = facnew/dot_product(p,ap)
call daxpy(nefc,alpha,p,1,x,1)