From 1ef2a0088753bc135cf22d3cc3c5811fdcc16eef Mon Sep 17 00:00:00 2001 From: Edoardo Apra Date: Wed, 22 Apr 2015 21:10:44 +0000 Subject: [PATCH] brought back (as alternative) the fmm based esp that was replaced after nwchem 6.1 --- src/esp/GNUmakefile | 2 + src/esp/esp.F | 25 ++++-- src/esp/esp_common.fh | 3 +- src/esp/esp_esp.F | 3 +- src/esp/esp_fastesp.F | 188 ++++++++++++++++++++++++++++++++++++++++++ src/esp/esp_rtdbin.F | 2 + 6 files changed, 214 insertions(+), 9 deletions(-) create mode 100644 src/esp/esp_fastesp.F diff --git a/src/esp/GNUmakefile b/src/esp/GNUmakefile index dc7bf0a2d4..cebee68286 100644 --- a/src/esp/GNUmakefile +++ b/src/esp/GNUmakefile @@ -7,6 +7,8 @@ esp_denmat.o\ esp_elemom.o\ esp_esp.o\ + esp_fastesp.o\ + esp_getden.o\ esp_fit.o\ esp_getg.o\ esp_getx.o\ diff --git a/src/esp/esp.F b/src/esp/esp.F index 6866714d0f..65fcb3e6a7 100644 --- a/src/esp/esp.F +++ b/src/esp/esp.F @@ -66,13 +66,24 @@ c c calculate the electrostatic potential on the grid c if(newgrd.ne.0.or.denmat.ne.0) then - call esp_esp(dbl_mb(ieptr), !ene - + dbl_mb(ioptr), !occ - + dbl_mb(idptr), !dns - + dbl_mb(ixptr), !x - + dbl_mb(iqptr), !q - + dbl_mb(igptr), !grid - + dbl_mb(ivptr)) !val + write(6,*) 'FASTESP ',fastesp + if(fastesp) then + call esp_fastesp(dbl_mb(ieptr), !ene + + dbl_mb(ioptr), !occ + + dbl_mb(idptr), !dns + + dbl_mb(ixptr), !x + + dbl_mb(iqptr), !q + + dbl_mb(igptr), !grid + + dbl_mb(ivptr)) !val + else + call esp_esp(dbl_mb(ieptr), !ene + + dbl_mb(ioptr), !occ + + dbl_mb(idptr), !dns + + dbl_mb(ixptr), !x + + dbl_mb(iqptr), !q + + dbl_mb(igptr), !grid + + dbl_mb(ivptr)) !val + endif end if ! newgrd, denmat c c do the fitting diff --git a/src/esp/esp_common.fh b/src/esp/esp_common.fh index 1310573575..adbfdee6e9 100644 --- a/src/esp/esp_common.fh +++ b/src/esp/esp_common.fh @@ -29,13 +29,14 @@ c integer ncon,icon(2500),irest,niter,maxit,ihfree,ngto,ngfr,idump integer lfngrd,nf,mbf,lscrn,newgrd,denmat,lfner,lfnplt,lfnxyz, + lfnchg + logical fastesp real*8 factor,probe,scrtol,qh common/esp_atm/rcut,spac,xmin,xmax,charge,chgnuc,chgele,radius, + rcon,resa,resb,resc,resd,toler,dipmom,rms,quamom,factor,probe, + scrtol,qh, + natoms,nbf,mgrid,mxgrid,ngrid,nset,nmo,nocc,nconst,ndim, + ncon,icon,irest,niter,maxit,ihfree,ngto,ngfr,idump,lfngrd,nf,mbf, - + lscrn,newgrd,denmat,lfner,lfnplt,lfnxyz,lfnchg + + lscrn,newgrd,denmat,lfner,lfnplt,lfnxyz,lfnchg,fastesp c c esp_chr : character variables c diff --git a/src/esp/esp_esp.F b/src/esp/esp_esp.F index 21d8f9e186..5c51ede188 100644 --- a/src/esp/esp_esp.F +++ b/src/esp/esp_esp.F @@ -51,7 +51,8 @@ c estimate mem. req. c call int_init_1eelec(maxbuf,maxscr,ibasis,0,ngrid) avail=ma_inquire_avail(mt_dbl)*8/10 - no_batches=nint(dble(maxbuf+maxscr)/dble(avail))+1 + no_batches=nint(dble(maxbuf+maxscr)/dble(avail)) + no_batches=max(no_batches,1) ngrid_batches=ngrid/no_batches #ifdef DEBUG if(me.eq.0) then diff --git a/src/esp/esp_fastesp.F b/src/esp/esp_fastesp.F new file mode 100644 index 0000000000..37e4cbafbd --- /dev/null +++ b/src/esp/esp_fastesp.F @@ -0,0 +1,188 @@ + subroutine esp_fastesp(ene,occ,dns,x,q,grid,val) +c +c $Id: esp_esp.F,v 1.22 2003-10-20 21:50:48 d3j191 Exp $ +c + implicit none +c +#include "esp_params.fh" +#include "esp_common.fh" +#include "global.fh" +#include "bas.fh" +#include "rtdb.fh" +#include "geom.fh" +#include "mafdecls.fh" +#include "msgids.fh" +#include "util.fh" +#include "schwarz.fh" +c + integer ga_create_atom_blocked + external ga_create_atom_blocked +c + real*8 ene(2*nbf),occ(2*nbf),dns(mbf,mbf) + real*8 x(3,natoms),q(natoms) + real*8 grid(3,mxgrid),val(mxgrid) +c + integer l_bt,ibtptr,l_st,istptr +c + integer i,j,ish,jsh,ix + integer ilo,ihi,jlo,jhi + real*8 dist,vt + integer imin,nval,jshi + real*8 dmin,fact + character*10 today,now + logical ldummy(1) +c +c integer lmax2, lencart, lenprod +c parameter (lmax2 = 2*3) ! Twice max angular momentum +c parameter (lencart = ((lmax2+1)*(lmax2+2))/2) +c parameter (lenprod = ((lmax2+1)*(lmax2+2)*(lmax2+3))/6) +c real*8 dwk(lencart,-lmax2:lmax2,0:lmax2), +c $ diwk(lenprod, -lmax2:lmax2,0:lmax2) +c + real*8 dwk(2548),diwk(7644) + integer npts + if(bas_is_spherical(ibasis)) + + call md_abort('fastesp does not work with spherical basis ',me) +c + call xlm_init() + call anl_fit_init() + call xlm_coeff_inv(6,dwk,diwk) +c +c calculate density matrix, occupation and energies +c + call esp_denmat(occ,ene) +c +c allocate memory for the integrals +c + if(.not.ma_push_get(mt_dbl,max1e,'bt',l_bt,ibtptr)) + + call md_abort('esp: ma_push_get b failed',me) + if(.not.ma_push_get(mt_dbl,mem1 ,'st',l_st,istptr)) + + call md_abort('esp: ma_push_get s failed',me) +c + + if(lscrn.ne.0) call schwarz_init(igeom,ibasis) +c +c get electrostatic potential on the grid points +c + do 2 i=1,ngrid + val(i)=0.0d0 + 2 continue +c + jshi=0 + if(me.eq.0) then + call util_file_name('er',.false.,.false.,erfil) + open(unit=lfner,file=erfil,form='unformatted',status='old', + + err=12) + read(lfner) jshi,nval + if(nval.ne.ngrid) call md_abort('Incompatible er file',0) + read(lfner) (val(i),i=1,ngrid) + close(unit=lfner) + 12 continue + write(lfnout,1000) jshi + 1000 format(/,' Recovering from shell ',i6,/) + endif + call ga_brdcst(mre_005,jshi,ma_sizeof(mt_int,1,mt_byte),0) +c + do 3 jsh=jshi+1,nshell + if(.not.bas_cn2bfr(ibasis,jsh,jlo,jhi)) + + call md_abort('esp: bas_cn2bfr failed',me) + do 4 ish=1,jsh + if(lscrn.ne.0) then + if(schwarz_shell(ish,jsh).lt.scrtol) goto 4 + endif + fact=1.0d0 + if(ish.ne.jsh) fact=2.0d0 + if(.not.bas_cn2bfr(ibasis,ish,ilo,ihi)) + + call md_abort('esp: bas_cn2bfr failed',me) + npts=ngto-ngfr+1 + call esp_getden(dns,jlo,jhi,ilo,ihi) + call potential_from_shell_pair(ibasis,jsh,ish,dns,mbf, + + scrtol,dwk,diwk,6,npts,grid(1,ngfr),val(ngfr),fact, + + .false.,ldummy,1) + 4 continue + call ga_dgop(mre_003,val,ngrid,'+') + if(me.eq.0) then + open(unit=lfner,file=erfil,form='unformatted',status='unknown') + rewind(unit=lfner) + write(lfner) jsh,ngrid + write(lfner) (val(i),i=1,ngrid) + close(unit=lfner) + if(util_print('output',print_medium)) then + call swatch(today,now) + write(lfnout,1001) today,now,jsh,nshell + 1001 format(2a10,' er shell ',i6,' of ',i6) + call util_flush(lfnout) + endif + else + do 13 i=1,ngrid + val(i)=zero + 13 continue + endif + 3 continue +c + if(me.eq.0) then + open(unit=lfner,file=erfil,form='unformatted',status='unknown') + close(unit=lfner,status='delete') + write(lfnout,1002) + 1002 format(/,' Recovery file deleted',/) + endif +c + do 8 i=ngfr,ngto + vt=val(i) + imin=0 + do 9 ix=1,natoms + dist=sqrt((grid(1,i)-x(1,ix))*(grid(1,i)-x(1,ix))+ + + (grid(2,i)-x(2,ix))*(grid(2,i)-x(2,ix))+ + + (grid(3,i)-x(3,ix))*(grid(3,i)-x(3,ix))) + if(imin.eq.0) then + dmin=dist + imin=ix + endif + if(dmin.gt.dist) then + dmin=dist + imin=ix + endif + val(i)=val(i)+q(ix)/dist + 9 continue + if(me.eq.0.and.util_print('where',print_debug)) then + if(i.eq.1) write(lfnout,1003) + 1003 format(//,' Electrostatic potential on grid',/) + write(lfnout,1004) imin,i,(grid(j,i),j=1,3),val(i)-vt,vt,val(i) + 1004 format(2i5,3f8.4,3f12.6) + endif + 8 continue +c + if(lscrn.ne.0) call schwarz_tidy() +c +c global sum of value array +c + if(np.gt.1) then + call ga_dgop(mre_006,val,ngrid,'+') + endif +c +c deallocate memory for the integrals +c + if(.not.ga_destroy(lg_d)) + + call md_abort('esp_denmat: ga_destroy lg_d failed',me) + if(.not.ma_pop_stack(l_st)) + + call md_abort('esp: ma_pop_stack st failed',me) + if(.not.ma_pop_stack(l_bt)) + + call md_abort('esp: ma_pop_stack bt failed',me) +c +c integral termination +c + call int_terminate() +c + if(me.eq.0) then + call util_file_name('grid',.false.,.false.,grdfil) + open(unit=lfngrd,file=grdfil,form='formatted',status='unknown') + rewind(lfngrd) + write(lfngrd,'(i10,f20.10)') ngrid,charge + do 10 i=1,ngrid + write(lfngrd,'(4f20.10)') (grid(j,i),j=1,3),val(i) + 10 continue + close(unit=lfngrd) + endif +c + return + end diff --git a/src/esp/esp_rtdbin.F b/src/esp/esp_rtdbin.F index 3f70ce625f..34206b2494 100644 --- a/src/esp/esp_rtdbin.F +++ b/src/esp/esp_rtdbin.F @@ -190,6 +190,8 @@ c if(.not.rtdb_get(irtdb,'esp:qh',mt_dbl,1,qh)) qh=0.0d0 if(.not.rtdb_get(irtdb,'esp:newgrd',mt_int,1,newgrd)) newgrd=0 if(.not.rtdb_get(irtdb,'esp:denmat',mt_int,1,denmat)) denmat=0 + if(.not.rtdb_get(irtdb,'esp:fastest',mt_log,1,fastesp)) + C fastesp=.false. c c convert rcut and spac from nm to au c