1 subroutine hnd_elpmap(rtdb,basis,geom) 2c 3c $Id$ 4c 5c This routine calculates the electrostatic potential 6c for a given density at the atomic positions. 7c 8 implicit none 9#include "errquit.fh" 10#include "global.fh" 11#include "mafdecls.fh" 12#include "nwc_const.fh" 13#include "stdio.fh" 14#include "geom.fh" 15#include "rtdb.fh" 16#include "cosmo.fh" 17#include "bas.fh" 18#include "msgids.fh" 19c 20 integer rtdb ! [Input] rtdb 21 integer basis ! [Input] Basis set 22 integer geom ! [Input] Geometry 23c 24 character*2 symbol 25 character*16 element, at_tag 26 integer iat, atn, nat, i 27 integer l_xyzpt, k_xyzpt, l_zanpt, k_zanpt, l_epot, k_epot 28 integer nefc, l_efcc, k_efcc, l_efcz, k_efcz 29 integer g_dens(3),ndens,nclosed(2),nopen(2),nvirt(2) 30 character*3 scftyp 31 double precision xp, yp, zp, xn, yn, zn, zan 32 double precision elpotn 33 double precision rr 34 double precision enuc 35c bq variables (MV) 36 logical dobq 37 integer bq_ncent 38 integer i_cbq 39 integer i_qbq 40 double precision elpotbq 41c property grid points variables (MV) 42 integer h_prp_c,i_prp_c 43 integer ma_prp_type 44 integer nprp 45 character*26 prp_date 46 logical do_points 47 logical do_output 48 integer l_tepot,k_tepot 49c 50 character*30 theory 51 integer nbf 52 integer ga_create_atom_blocked 53 external ga_create_atom_blocked 54 logical ao_1prdm_read 55 external ao_1prdm_read 56 logical ocube 57 integer nsp(3) 58 integer un_grid,avail_ma 59 integer npasses,i_pass 60 integer iptr,nprp_pass,ngrid(3) 61 character*255 cube_name 62 character*255 dmat_file 63c 64c 65c Initialize integrals 66c 67 call int_init(rtdb,1, basis) 68 call schwarz_init(geom, basis) 69c 70c Get density matrix 71c 72 if(.not.rtdb_cget(rtdb,'task:theory',1,theory)) 73 + call errquit('task: no task input for theory?',0, RTDB_ERR) 74 if (.not. rtdb_cget(rtdb,'scf:scftype',1,scftyp)) scftyp = "RHF" 75 ndens = 1 76 if (scftyp.eq.'UHF') ndens = 3 77c 78c if (ga_nodeid().eq.0) write(luout,*) "scftyp:",scftyp 79c if (ga_nodeid().eq.0) write(luout,*) "theory:",theory 80c 81c Read density matrix from a file 82c 83 if (theory.eq.'tce'.or.theory.eq.'tddft') then 84c 85 do i = 1, ndens 86 g_dens(i) = ga_create_atom_blocked(geom,basis,'density matrix') 87 call ga_zero(g_dens(i)) 88 enddo 89c 90 if (.not. bas_numbf(basis,nbf)) 91 & call errquit('hnd_elfmap: could not get nbf',0, BASIS_ERR) 92c 93 call util_file_name('dmat',.false.,.false.,dmat_file) ! get filename 94 if (ga_nodeid().eq.0) 95 & write(luout,*) "Reading file from:",dmat_file 96 if(.not.ao_1prdm_read(nbf,g_dens(ndens),dmat_file)) 97 & call errquit('hnd_elfmap: ao_1prdm_read failed',0,0) 98c 99 else ! generate from movecs for scf or dft 100 call hnd_prp_get_dens(rtdb,geom,basis,g_dens,ndens,scftyp, 101 & nclosed,nopen,nvirt) 102 endif 103c 104 call ecce_print_module_entry('Elpoten') 105c 106c ----- define points for calculation ----- 107c 1. grid points (not active) 108c 2. nuclei 109c 3. center of mass (not active) 110c 111 if (.not.geom_ncent(geom,nat)) call 112 & errquit('hnd_elpmap: geom_ncent',911,GEOM_ERR) 113c 114 if (.not. ma_push_get(mt_dbl,3*nat,'xyz pnt',l_xyzpt,k_xyzpt)) 115 & call errquit('hnd_elpmap: ma failed',1,MA_ERR) 116 if (.not. ma_push_get(mt_dbl,nat,'zan pnt',l_zanpt,k_zanpt)) 117 & call errquit('hnd_elpmap: ma failed',2,MA_ERR) 118c 119 do 30 iat=1,nat 120 if(.not.geom_cent_get(geom,iat,at_tag,dbl_mb(k_xyzpt+3*(iat-1)), 121 & dbl_mb(k_zanpt+iat-1))) call 122 & errquit('hnd_elpmap: geom_cent_get',911,GEOM_ERR) 123 30 continue 124c 125 if(.not.rtdb_get(rtdb,'prop:cubefile',mt_log,1,ocube)) 126 + ocube = .false. 127 128 if(ocube) then 129 if (.not. rtdb_cget(rtdb, "prop:grid:output",1,cube_name)) 130 > call util_file_prefix("elp.cube",cube_name) 131 end if 132cc 133c define points for the calculation now 134c either custom grid or (default) nuclei positions (M.V.) 135c ------------------------------------------------- 136 if(ocube) then 137 call prop_grid_initialize(rtdb,nat,dbl_mb(k_xyzpt)) 138 call prop_grid_get_r_ptr(nprp,i_prp_c,ngrid) 139 do_points = .false. 140 do_output = .false. 141 else if(rtdb_get_info(rtdb, "prop:xyz", ma_prp_type, 142 > nprp, prp_date)) then 143 nprp = nprp/3 144 if (.not. ma_push_get(mt_dbl,3*nprp,'prop:xyz',h_prp_c,i_prp_c)) 145 & call errquit('hnd_elpmap: prop:xyz',911,MA_ERR) 146 if (.not. rtdb_get(rtdb,'prop:xyz',mt_dbl, 147 > 3*nprp,dbl_mb(i_prp_c))) 148 & call errquit('hnd_elpmap: prop:xyz failed',911,RTDB_ERR) 149 do_points = .true. 150 do_output = .true. 151 else 152 nprp = nat 153 if (.not. ma_push_get(mt_dbl,3*nat,'prop:xyz',h_prp_c,i_prp_c)) 154 & call errquit('hnd_elpmap: ma failed',3,MA_ERR) 155 call dcopy(3*nat,dbl_mb(k_xyzpt),1,dbl_mb(i_prp_c),1) 156 do_points = .false. 157 do_output = .true. 158 end if 159c 160#ifndef OLDCUBE 161 if(ocube) 162 c call prop_grid_writecubehead(geom,un_grid,cube_name) 163#endif 164c 165c check if enough mem is available 166c 167 avail_ma = ma_inquire_avail(mt_dbl)*0.9d0 168 npasses=nprp*2d0/avail_ma + 1d0 169c 170 call ga_igop(msg_efgs_col-4,npasses,1,'max') 171#ifdef OLDCUBE 172 if(npasses.gt.1) call errquit('hnd_elpmap: ma failed',0,MA_ERR) 173#endif 174 if(npasses.gt.1.and.do_points) 175 C call errquit('hnd_elpmap: ma failed',0,MA_ERR) 176 nprp_pass=nprp/npasses 177 if(npasses.gt.1) then 178c 179c enforce nprp_pass to be a multiple of ngrid(3) to keep cube happy 180c 181 nprp_pass=(nprp_pass/ngrid(3))*ngrid(3) 182 npasses=nprp/nprp_pass 183 if(ga_nodeid().eq.0) then 184 write(6,'(a,i5,a,i10,a,i19)') 185 H ' hnd_elpmap: npasses:',npasses, 186 W ' nprp_pass:',nprp_pass,' nprp:',nprp 187 endif 188 endif 189 iptr=0 190 do i_pass=1,npasses 191 iptr=(i_pass-1)*nprp_pass 192 if(i_pass.eq.npasses) nprp_pass=nprp-iptr 193 194 if (.not. ma_push_get(mt_dbl,nprp_pass,'epot pnt',l_epot,k_epot)) 195 & call errquit('hnd_elpmap: ma failed',4,MA_ERR) 196c 197c total electric field array (M.V.) 198c -------------------------------- 199 if (.not. ma_push_get(mt_dbl,nprp_pass,'totepot',l_tepot,k_tepot)) 200 & call errquit('hnd_elpmap: ma failed',5,MA_ERR) 201c 202cc 203c ----- calculate electronic contribution at all points ----- 204c 205 call hnd_elfcon(basis,geom,g_dens(ndens),dbl_mb(i_prp_c+iptr*3), 206 N nprp_pass, dbl_mb(k_epot),0) 207c 208 dobq = .false. 209 bq_ncent = 0 210 i_cbq = 0 211 i_qbq = 0 212 if(geom_extbq_on()) then 213 dobq = .true. 214 bq_ncent = geom_extbq_ncenter() 215 i_cbq = geom_extbq_coord() 216 i_qbq = geom_extbq_charge() 217 end if 218 219c 220c ----- collect and output results of all points ----- 221c 222 if (ga_nodeid().gt.0) goto 300 223c 224c ----- calculate electrostatic potential ----- 225c 226 if(do_output) then 227 if (ga_nodeid().eq.0) write(luout,9999) 228 if (ga_nodeid().eq.0) write(luout,9994) 229c 230c 231 if(dobq) then 232 write(luout,9992) 233 else 234 write(luout,9997) 235 end if 236 end if 237 do 230 iat=1,nprp_pass 238 xp = dbl_mb(iptr*3+i_prp_c +3*(iat-1)) 239 yp = dbl_mb(iptr*3+i_prp_c+1+3*(iat-1)) 240 zp = dbl_mb(iptr*3+i_prp_c+2+3*(iat-1)) 241c 242c ----- add nuclear contribution ----- 243c 244 elpotn = -dbl_mb(k_epot+iat-1) 245 elpotbq = 0.0d0 246 enuc = 0.0 247 do 210 i = 1,nat 248 xn = dbl_mb(k_xyzpt +3*(i-1)) - xp 249 yn = dbl_mb(k_xyzpt+1+3*(i-1)) - yp 250 zn = dbl_mb(k_xyzpt+2+3*(i-1)) - zp 251 zan = dbl_mb(k_zanpt+i-1) 252 rr = sqrt(xn*xn + yn*yn + zn*zn) 253 if (rr.lt.1.0d-3) go to 210 254 elpotn = elpotn + zan/rr 255 enuc = enuc + zan/rr 256 210 continue 257c 258c ----- form -efc- contribution ----- 259c from cosmo point charges !!!! 260c 261 if (cosmo_last) then 262 if (.not.rtdb_get(rtdb,'cosmo:nefc',mt_int,1,nefc)) 263 & call errquit('hnd_elpmap: rtdb get failed for nefc ',911, 264 & RTDB_ERR) 265 if (.not.ma_push_get(mt_dbl,nefc*3,'efcc',l_efcc,k_efcc)) 266 & call errquit('hnd_elpmap: malloc k_efcc fail',911,ma_err) 267 if (.not.ma_push_get(mt_dbl,nefc,'efcz',l_efcz,k_efcz)) 268 & call errquit('hnd_elpmap: malloc k_efcz fail',911,ma_err) 269 if (.not.rtdb_get(rtdb,'cosmo:efcc',mt_dbl,3*nefc, 270 & dbl_mb(k_efcc))) call 271 & errquit('hnd_elpmap: rtdb get failed efcc',912,rtdb_err) 272 if (.not.rtdb_get(rtdb,'cosmo:efcz',mt_dbl,nefc, 273 & dbl_mb(k_efcz))) call 274 & errquit('hnd_elpmap: rtdb get failed efcz',913,rtdb_err) 275 do 220 i = 1,nefc 276 xn = dbl_mb(k_efcc+3*(i-1) ) - xp 277 yn = dbl_mb(k_efcc+3*(i-1)+1) - yp 278 zn = dbl_mb(k_efcc+3*(i-1)+2) - zp 279 rr = sqrt(xn*xn + yn*yn + zn*zn) 280 if (rr.lt.1.0d-3) then 281 write(luout,9993) xp,yp,zp,i 282 go to 220 283 endif 284 elpotn = elpotn + dbl_mb(k_efcz+i-1)/rr 285 220 continue 286 if (.not.ma_chop_stack(l_efcc)) call 287 & errquit('hnd_elpmap: chop stack l_efcc',913,ma_err) 288 endif 289c adding external bq contributions(MV) 290c ---------------------------------- 291 if (dobq) then 292 do i = 1,bq_ncent 293 xn = dbl_mb(i_cbq+3*(i-1) ) - xp 294 yn = dbl_mb(i_cbq+3*(i-1)+1) - yp 295 zn = dbl_mb(i_cbq+3*(i-1)+2) - zp 296 rr = sqrt(xn*xn + yn*yn + zn*zn) 297 elpotbq = elpotbq+dbl_mb(i_qbq+i-1)/rr 298 end do 299 end if 300 elpotn = elpotn + elpotbq 301 dbl_mb(k_tepot+iat-1) = elpotn 302c end of external bq contributions 303c ------------------------------- 304 if(do_output) then 305 if(.not.do_points) then 306 if (.not. geom_cent_tag(geom,iat,at_tag)) call 307 & errquit('hnd_elpmap: geom_cent_tag failed',0,GEOM_ERR) 308c geom_tag_to_element returns false for Bq elements (MV) 309c ----------------------------------------------------- 310 if (.not. geom_tag_to_element(at_tag,symbol,element,atn)) then 311 if(symbol.ne."bq") call 312 & errquit('hnd_elpmap: geom_tag_to_element failed',0,GEOM_ERR) 313 end if 314 end if 315c 316 if(do_points) then 317 symbol = "pt" 318 if(dobq) then 319 write(luout,9991) iat,symbol,xp,yp,zp,elpotn, 320 & dbl_mb(k_epot+iat-1),elpotbq 321 else 322 write(luout,9995) iat,symbol,xp,yp,zp,enuc, 323 & dbl_mb(k_epot+iat-1) 324 end if 325 else 326 if(dobq) then 327 write(luout,9991) iat,symbol,xp,yp,zp,elpotn, 328 & dbl_mb(k_epot+iat-1),elpotbq 329 else 330 write(luout,9995) iat,symbol,xp,yp,zp,elpotn, 331 & dbl_mb(k_epot+iat-1) 332 end if 333 end if 334 end if 335c 336c ----- store ecce data ----- 337c 338 call ecce_print1_char('atom name',symbol,1) 339 call ecce_print1('Electrostatic pot',MT_DBL,elpotn,1) 340 call ecce_print1('Diamagnetic shield',MT_DBL, 341 & dbl_mb(k_epot+iat-1),1) 342c 343 230 continue ! Assembling and printing next atom 344c 345 call ecce_print_module_exit('Elpoten','ok') 346 call util_flush(luout) 347c 348c ----- release memory block ----- 349c 350 300 call ga_sync() 351 352 if(ga_nodeid().eq.0.and.ocube) then 353 call util_print_centered(luout, 354 > "writing esp potential to "//cube_name,.true.,.true.) 355#ifdef OLDCUBE 356 call prop_grid_write_cube(geom,nprp,dbl_mb(k_tepot),cube_name) 357#else 358c call prop_grid_writecubehead(geom,un_grid,cube_name) 359 call prop_grid_writecubegrid(nprp_pass,dbl_mb(k_tepot), 360 F un_grid) 361#endif 362 end if 363c 364c if custom grid is requested save final electric 365c into rtdb(M.V.) 366c ----------------------------------------------- 367 if(do_points) then 368 if (.not. rtdb_put(rtdb,'prop:epot_xyz',mt_dbl, 369 > nprp,dbl_mb(k_tepot))) 370 & call errquit('hnd_elpmap: epot_xyz failed',0,RTDB_ERR) 371 end if 372cc 373c ------- Deallocate MA memory ------ 374c 375 if (.not.ma_pop_stack(l_tepot)) call errquit 376 & ('hnd_elpmap, ma_pop_stack of l_tepot failed',911,MA_ERR) 377 if (.not.ma_pop_stack(l_epot)) call errquit 378 & ('hnd_elpmap, ma_pop_stack of l_epot failed',911,MA_ERR) 379 enddo ! i_pass 380 if(.not.ocube) then 381 if (.not.ma_pop_stack(h_prp_c)) call errquit 382 & ('hnd_elpmap, ma_pop_stack of h_prp_c failed',911,MA_ERR) 383 else 384 close(un_grid) 385 call prop_grid_destroy() 386 endif 387 if (.not.ma_pop_stack(l_zanpt)) call errquit 388 & ('hnd_elpmap, ma_pop_stack of l_zanpt failed',911,MA_ERR) 389 if (.not.ma_pop_stack(l_xyzpt)) call errquit 390 & ('hnd_elpmap, ma_pop_stack of l_xyzpt failed',911,MA_ERR) 391c 392 do i = 1, ndens 393 if (.not.ga_destroy(g_dens(i))) call 394 & errquit('elpmap: ga_destroy failed g_dens',0,GA_ERR) 395 enddo 396c 397c Terminate integrals 398c 399 call schwarz_tidy() 400 call int_terminate() 401c 402 return 403 9999 format(/,10x,45(1H-), 404 1 /,10x,'Electrostatic potential/diamagnetic shielding', 405 2 /,10x,45(1H-),/) 406 9998 format(' Not enough core in -elpmap-') 407 9997 format(3x,'Point',6x,'X',9x,'Y',9x,'Z',5x,'Total Potential(a.u.)', 408 1 3x,'Diamagnetic shielding(a.u.)') 409 9996 format(' --- Warning - electrostatic potential at ', 410 1 3f10.5,' . contribution from nucleus ',i3,' ignored') 411 9995 format(i5,1x,a2,3F10.5,f15.6,6x,f15.6) 412 9994 format(' 1 a.u. = 9.07618 esu/cm ( or statvolts ) ') 413 9993 format(' --- Warning - electrostatic potential at ', 414 1 3f10.5,' . contribution from -efc- ',i3,' ignored') 415c 416 9992 format(3x,'Point',6x,'X',9x,'Y',9x,'Z',5x,'Total Potential(a.u.)', 417 1 3x,'Diamagnetic shielding(a.u.)', 418 2 3x,'External Bq potential') 419 9991 format(i5,1x,a2,3F10.5,f15.6,6x,f15.6,12x,f15.6) 420 421 END 422