diff --git a/src/nwdft/grid/grid_checks.F b/src/nwdft/grid/grid_checks.F
index a97269f8a1..5fce7a1d9f 100644
--- a/src/nwdft/grid/grid_checks.F
+++ b/src/nwdft/grid/grid_checks.F
@@ -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
diff --git a/src/nwdft/grid/grid_setang.F b/src/nwdft/grid/grid_setang.F
index 4b959477cf..9670e051b2 100644
--- a/src/nwdft/grid/grid_setang.F
+++ b/src/nwdft/grid/grid_setang.F
@@ -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
diff --git a/src/nwdft/xc/xc_lebedev.F b/src/nwdft/xc/xc_lebedev.F
index 67200cb038..d39aa2f853 100644
--- a/src/nwdft/xc/xc_lebedev.F
+++ b/src/nwdft/xc/xc_lebedev.F
@@ -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)
diff --git a/src/solvation/GNUmakefile b/src/solvation/GNUmakefile
index efc8a32413..118cf4f349 100644
--- a/src/solvation/GNUmakefile
+++ b/src/solvation/GNUmakefile
@@ -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
diff --git a/src/solvation/cosmo_charges.F b/src/solvation/cosmo_charges.F
index e578b85c30..ab19c4130a 100644
--- a/src/solvation/cosmo_charges.F
+++ b/src/solvation/cosmo_charges.F
@@ -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 -----
diff --git a/src/solvation/cosmo_initialize.F b/src/solvation/cosmo_initialize.F
index ab2e7c030e..d9d95465c9 100644
--- a/src/solvation/cosmo_initialize.F
+++ b/src/solvation/cosmo_initialize.F
@@ -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> @}
diff --git a/src/solvation/cosmo_input.F b/src/solvation/cosmo_input.F
index c88b2f40da..96d5127113 100644
--- a/src/solvation/cosmo_input.F
+++ b/src/solvation/cosmo_input.F
@@ -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
diff --git a/src/solvation/cosmo_lebedev.F b/src/solvation/cosmo_lebedev.F
new file mode 100644
index 0000000000..757889faf4
--- /dev/null
+++ b/src/solvation/cosmo_lebedev.F
@@ -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
diff --git a/src/solvation/cosmo_params.fh b/src/solvation/cosmo_params.fh
index debd050386..c47d76ee61 100644
--- a/src/solvation/cosmo_params.fh
+++ b/src/solvation/cosmo_params.fh
@@ -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
diff --git a/src/solvation/grad_hnd_cos.F b/src/solvation/grad_hnd_cos.F
index 6a93892abd..10c7e5cd6c 100644
--- a/src/solvation/grad_hnd_cos.F
+++ b/src/solvation/grad_hnd_cos.F
@@ -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)
diff --git a/src/solvation/hnd_cosmo_lib.F b/src/solvation/hnd_cosmo_lib.F
index 3e98a0c681..81886df7a9 100644
--- a/src/solvation/hnd_cosmo_lib.F
+++ b/src/solvation/hnd_cosmo_lib.F
@@ -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>
C> 10.1111/1468-0262.00296.
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>
C> 10.1021/jp992097l.
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> J. Chem. Soc., Perkin Trans. 2, 1993, pp 799-805, DOI:
C>
C> 10.1039/P29930000799.
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)