1 subroutine hnd_giaox_zora(rtdb,basis,geom) 2c $Id$ 3c 4 implicit none 5c 6#include "errquit.fh" 7#include "global.fh" 8#include "mafdecls.fh" 9#include "msgids.fh" 10#include "geom.fh" 11#include "rtdb.fh" 12#include "bas.fh" 13#include "stdio.fh" 14#include "apiP.fh" 15#include "prop.fh" 16#include "bgj.fh" 17#include "case.fh" 18#include "inp.fh" 19#include "zora.fh" 20c 21 integer rtdb ! [input] rtdb handle 22 integer basis ! [input] basis handle 23 integer geom ! [input] geometry handle 24 integer nclosed(2), nopen(2), nvirt(2), ndens, nbf, nmo 25 integer ixy, ix, iy, iatom, iocc, ifld, ioff 26 integer alo(3), ahi(3), blo(3), bhi(3), clo(3), chi(3) 27 integer dlo(3), dhi(3) 28 integer l_occ, k_occ, l_eval, k_eval 29 integer l_dia, k_dia, l_para, k_para 30 integer l_xyz, k_xyz, l_zan, k_zan 31 integer l_AtNr, k_AtNr,type_NMR 32 integer l_fukui,k_fukui ! to include Fukui term 33 ! allows dia,para tensors 34 ! to be origin invariant separately 35 integer g_dens(3), g_s10, g_h01, g_h11, 36 & g_d1, g_rhs, g_rhs0,g_fock, g_u 37 integer g_fock_Coul,g_fock_Exch, 38 & g_rhs_Coul,g_rhs_Exch,g_rhs_noJK, 39 & g_rhs_eSji,g_rhs_1e 40 integer g_d1_oo,g_d1_ov ! split occ-occ,occ-virt contrib to 41 ! perturbed density matrix, P10 01-26-11 42 integer g_d1_ov_Coul,g_d1_ov_Exch,g_d1_ov_noJK, 43 & g_d1_ov_1e,g_d1_ov_eSji,g_d1_ov_ExchSFB 44 45 integer vectors(2), geomnew, i, j, ij, g_xc(3) 46 integer vectors_scl(2) ! to be scaled with g_CiFull 47 double precision atn, tol2e, val, isotr, aniso 48 double precision val_oo,val_ov 49 double precision a(6),axs(3,3),eig(3),xfac 50 double precision jfac(12),kfac(12) 51 character*255 zorafilename 52 character*3 scftyp 53 character*16 tag 54 character*32 element 55 character*256 cphf_rhs, cphf_sol 56 character*2 symbol 57 integer ld(2),cbuf ! FA-08-19-10 58 logical cphf2, file_write_ga, file_read_ga, cphf 59 external cphf2, file_write_ga, file_read_ga, cphf 60 external giao_aotomo,hnd_giasym,mat_transpose, ! FA-09-16-10 61 & dft_zoraNMR_read,get_slctd_atoms 62 logical dft_zoraNMR_read 63 double precision val1 64 logical oskel, status,debug_cs 65 double precision val_ov_Coul,val_ov_Exch,val_ov_noJK, 66 & val_ov_1e,val_ov_eSji 67 68 double precision ppm 69 data ppm /26.62566914d+00/ 70 data tol2e /1.0d-10/ 71 integer ic,npol 72c ----- Definitions for NLMO analysis ---- START 73 integer i1,j1,acc_vec,l_tvec,k_tvec,shldfile 74 integer g_tvec ! for nbo analysis 75c ----- Definitions for NLMO analysis ---- END 76 external get_prelim_fock,get_P10,add_H10,add_fock, 77 & get_P10_1,skip_cphf,get_par,get_dia, 78 & skip_cphf_JK,get_P10_JK,new_giao_2e_JK, 79 & get_par_JK, 80 & dft_zoraCPHF_write,dft_zoraCPHF_read 81 integer ntot,ispin,indA,indB,nocc(2),nind_jk, 82 & typeprop,nat_slc,nat 83 integer cdens,l_pararr,k_pararr,icalczora,type_nmrdata 84 integer ga_h01,npar_analysis 85 integer debug_giaox 86 integer shift,disp,nbf_dat,nat_dat 87 logical switch_nmrcs_analysis,switch_skip_cphf 88 integer g_AtNr1 89 integer ga_dia, ! OUTPUT 90 & ga_para1, ! OUTPUT 91 & ga_h01_num,! OUTPUT 92 & ga_Fji ! OUTPUT 93 logical skip_cphf_ev_shield 94 oskel=.false. 95c Note.- switch_nmrcs_analysis=.true. has to go together with 96c switch_skip_cphf=.true. 97c ----------------------------------- 98 debug_giaox=0 ! =1 for debugging print outs of matrices 99c debug_giaox=1 ! =1 for debugging print outs of matrices 100c switch_skip_cphf=.true. ! For skipping cphf or cpks 101 switch_skip_cphf=.false. ! For NOT skipping cphf or cpks 102 switch_nmrcs_analysis=.false. ! For default mode NO nmrcs analysis 103c switch_nmrcs_analysis=.true. ! For analysis of nmrcs tensor 104 if (ga_nodeid().eq.0.and.debug_giaox.eq.1) 105 & write(*,*) 'switch_skip_cphf=',switch_skip_cphf 106 if (ga_nodeid().eq.0.and.debug_giaox.eq.1) 107 & write(*,*) 'switch_nmrcs_analysis=',switch_nmrcs_analysis 108c 109c Print NMR shielding header 110c 111 if (ga_nodeid().eq.0) write(luout,9999) 112c 113c Current CPHF does not handle symmetry 114c Making C1 geometry and store it on rtdb 115c 116c If DFT get part of the exact exchange defined 117c 118 xfac = 1.0d0 119 if (use_theory.eq.'dft') xfac = bgj_kfac() 120 nind_jk=12 121 do ifld = 1,nind_jk 122 jfac(ifld) = 0.0d0 ! used in shell_fock_build() 123 kfac(ifld) = -1.0d0*xfac ! used in shell_fock_build() 124 if (ga_nodeid().eq.0.and.debug_giaox.eq.1) then 125 write(*,144) ifld,jfac(ifld),kfac(ifld) 126 144 format('(j,k)(',i3,')=(',f15.8,',',f15.8,')') 127 endif 128 enddo 129c 130c Integral initialization 131 call int_init(rtdb,1,basis) 132 call schwarz_init(geom,basis) 133 call hnd_giao_init(basis,1) 134 call scf_get_fock_param(rtdb,tol2e) 135c 136c Find out from rtdb which atoms we need to calculate shielding for 137c Get number of atoms (all or number from rtdb) 138c Get which atoms (all or some read from rtdb) 139c Allocate arrays which will hold atomic information (k_zan and k_xyz) 140 141 status = rtdb_parallel(.true.) 142c ------- Read (nat,atmnr) --------- START 143 status=geom_ncent(geom,nat) 144 if (.not.ma_alloc_get( 145 & mt_int,nat,'nmt tmp',l_AtNr,k_AtNr)) 146 & call errquit('hnd_giaox_zora: ma failed',0,MA_ERR) 147 typeprop=2 ! =1 EFG =2 Shieldings =3 Hyperfine 148 call get_slctd_atoms(nat_slc, ! out: selected atoms 149 & int_mb(k_AtNr), ! out: list of selected atom nr. 150 & nat, ! in : total nr atoms in molecule 151 & rtdb, ! in : rdt handle 152 & typeprop) ! in : =1,2,3=EFG,Shieldings,Hyperfine 153 if (ga_nodeid().eq.0.and.debug_giaox.eq.1) then 154 write(*,*) 'nat_slc=',nat_slc 155 do i=1,nat_slc 156 write(*,7) i,int_mb(k_AtNr+i-1) 157 7 format('atomnr(',i3,')=',i5) 158 enddo 159 endif 160c ------- Read (nat,atmnr) --------- END 161 if (.not. ma_push_get(mt_dbl,3*nat_slc,'nmr at',l_xyz,k_xyz)) 162 & call errquit('hnd_giaox_zora: ma_push_get failed k_xyz', 163 & 0,MA_ERR) 164 if (.not. ma_push_get(mt_dbl,nat_slc,'nmr zan',l_zan,k_zan)) 165 & call errquit('hnd_giaox_zora: ma_push_get failed k_zan', 166 & 0,MA_ERR) 167c 168c Try to read the atom list from rtdb. If it is not there, we still have the default list 169 do ixy = 0, nat_slc-1 170 if (.not. geom_cent_get(geom, int_mb(k_AtNr+ixy), tag, 171 & dbl_mb(k_xyz+3*ixy),dbl_mb(k_zan+ixy))) 172 & call errquit('hnd_giaox_zora: geom_cent_tag failed', 173 & 0, GEOM_ERR) 174 enddo 175c 176c Get Unperturbed MO vectors and eigenvalues 177c First allocate some memory for occupation numbers and eigenvalues 178 179 if (.not. bas_numbf(basis,nbf)) call 180 & errquit('hnd_giaox_zora: could not get nbf',0, BASIS_ERR) 181 182c ++++++ Reading shieldings data from file ++++++ START 183 if (do_zora) then ! ==========if-do_zora-START 184c Note.- lbl_nmrcs defined in zora.fh 185 call util_file_name(lbl_nmrcs,.false.,.false.,zorafilename) 186 icalczora = 0 ! initialize the flag 187 type_nmrdata=1 ! =1,2,3=shieldings,hyperfine,gshift 188 if (.not.dft_zoraNMR_read(zorafilename, 189 & type_nmrdata, 190 & nbf,nat_slc, 191 & g_AtNr1,ga_dia, 192 & ga_para1,ga_h01_num,ga_Fji)) icalczora=1 193c Note.- If I print the GAs here it gets freezed 194 endif ! ======================if-do_zora-END 195c ++++++ Reading shieldings data from file ++++++ END 196 if (debug_giaox.eq.1) then 197 if (ga_nodeid().eq.0) 198 & write(*,*) '----ga_para1---------- START' 199 call ga_print(ga_para1) 200 if (ga_nodeid().eq.0) 201 & write(*,*) '----ga_para1---------- END' 202 if (ga_nodeid().eq.0) 203 & write(*,*) '----ga_h01_num---------- START' 204 call ga_print(ga_h01_num) 205 if (ga_nodeid().eq.0) 206 & write(*,*) '----ga_h01_num---------- END' 207 if (ga_nodeid().eq.0) 208 & write(*,*) '----ga_Fji---------- START' 209 call ga_print(ga_Fji) 210 if (ga_nodeid().eq.0) 211 & write(*,*) '----ga_Fji---------- END' 212 endif 213 if (.not. ma_push_get(mt_dbl,2*nbf,'occ num',l_occ,k_occ)) call 214 & errquit('hnd_giaox_zora: ma_push_get failed k_occ',0,MA_ERR) 215 if (.not. ma_push_get(mt_dbl,2*nbf,'eigenval',l_eval,k_eval)) call 216 & errquit('hnd_giaox_zora: ma_push_get failed k_eval',0,MA_ERR) 217c 218 call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen, 219 & nvirt,scftyp,vectors,dbl_mb(k_occ), 220 & dbl_mb(k_eval),nmo) 221c 222c 223c ------ define npol ----- START 224 if (.not. rtdb_get(rtdb, 'dft:ipol', mt_int, 1, npol)) then 225 if (scftyp .eq. 'UHF') then 226 npol=2 227 else if (scftyp .eq. 'RHF') then 228 npol=1 229 endif 230 endif 231c ------ define npol ----- END 232c if (ga_nodeid().eq.0) 233c & write(*,*) 'npol=',npol 234 235 if (npol.eq.1) then 236 nocc(1)=nclosed(1) 237 nocc(2)=0 238 else if (npol.eq.2) then 239 nocc(1)=nopen(1) 240 nocc(2)=nopen(2) 241 endif 242 if (debug_giaox.eq.1) then 243 write(*,10) nocc(1),nocc(2), 244 & nopen(1),nopen(2), 245 & nclosed(1),nclosed(2), 246 & nvirt(1),nvirt(2),scftyp,nmo 247 10 format('nocc =(',i3,',',i3,') ', 248 & 'nopen=(',i3,',',i3,') ', 249 & 'nclos=(',i3,',',i3,') ', 250 & 'nvirt=(',i3,',',i3,') scftyp=',a,', nmo=',i3) 251 endif 252c ------ Store nopen in rtdb so that CPHF routine is happy ---- START 253c WARNING: For restricted calc nocc(1)=9 nocc(2)=0 <=== PROBLEM 254 if (npol.eq.2) then 255 if (.not. rtdb_put(rtdb, 'scf:nopen', 256 & MT_INT, 1, nocc(1)-nocc(2))) 257 * call errquit('hnd_giaox_zora:rtdbput nopen failed', 258 & nocc(1)-nocc(2), 259 & RTDB_ERR) 260 endif 261c ------ Store nopen in rtdb so that CPHF routine is happy ---- END 262c 263c Get Unperturbed Density Matrix 264 call hnd_prp_get_dens(rtdb,geom,basis, 265 & g_dens, ! out 266 & ndens,scftyp, 267 & nclosed,nopen,nvirt) 268c 269 ntot=0 270 do ispin=1,npol 271 ntot=ntot+nocc(ispin)*nvirt(ispin) 272 enddo 273 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs)) 274 & call errquit('hnd_giaox_zora: ga_create failed g_rhs',0,GA_ERR) 275 call ga_zero(g_rhs) 276 277 if (switch_nmrcs_analysis) then 278 if (ga_nodeid().eq.0) 279 & write(*,*) 'ENTER-define g_rhs(J,K,nJK,eSji' 280c ----- TEST: for testing Coulomb and Exchange contrib --- START 281 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_Coul)) 282 & call errquit('hnd_giaox_zora: ga_create failed g_rhsJ', 283 & 0,GA_ERR) 284 call ga_zero(g_rhs_Coul) 285 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_Exch)) 286 & call errquit('hnd_giaox_zora: ga_create failed g_rhsK', 287 & 0,GA_ERR) 288 call ga_zero(g_rhs_Exch) 289 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_noJK)) 290 & call errquit('hnd_giaox_zora: ga_create failed g_rhsnJK', 291 & 0,GA_ERR) 292 call ga_zero(g_rhs_noJK) 293 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_eSji)) 294 & call errquit('hnd_giaox_zora: ga_create failed g_rhseSji', 295 & 0,GA_ERR) 296 call ga_zero(g_rhs_eSji) 297 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs_1e)) 298 & call errquit('hnd_giaox_zora: ga_create failed g_rhs1e', 299 & 0,GA_ERR) 300 call ga_zero(g_rhs_1e) 301c ----- TEST: for testing Coulomb and Exchange contrib --- END 302 endif 303 304 if (debug_giaox.eq.1) then 305 write(*,*) 'after-creating g_rhs ...' 306 write(*,102) npol,nocc(1),nocc(2), 307 & nclosed(1),nclosed(2), 308 & nvirt(1),nvirt(2),scftyp,ntot 309 102 format('BEF pre-fock::npol=',i3,' nocc=(',i3,',',i3,') ', 310 & 'nclos=(',i3,',',i3,') ', 311 & 'nvirt=(',i3,',',i3,') scftyp=',a,' ntot=',i3) 312 endif 313 314 if (switch_nmrcs_analysis) then 315 call get_prelim_fock_debug( 316 & g_d1, ! out: 317 & g_rhs, ! out: rhs expression 318 & g_rhs0, ! out: to be used in get_d1() 319 & g_rhs_eSji, ! out: -Sji^k \epsilon_i 320 & vectors, ! in: MO coeffs 321 & dbl_mb(k_eval), ! in: energy vals 322 & dbl_mb(k_xyz), ! in: Nuclear positions (x,y,z) 323 & nat_slc, ! in: nr. selected nuclei (atoms) 324 & basis, ! in: basis handle 325 & nbf, ! in: nr. basis functions 326 & nmo, ! in: nr. MOs (occ+virt) 327 & npol, ! in: nr. of polarizations 328 & nocc, ! in: nr. occ MOs 329 & nvirt) ! in: nr. virtual MOs 330 else ! default 331 call get_prelim_fock(g_d1, ! out: 332 & g_rhs, ! out: rhs expression 333 & g_rhs0,! out: to be used in get_d1() 334 & vectors,dbl_mb(k_eval), 335 & dbl_mb(k_xyz),nat_slc,basis, 336 & nbf,nmo,npol,nocc,nvirt) 337 endif 338c 339 if (debug_giaox.eq.1) then 340 if (ga_nodeid().eq.0) 341 & write(*,*) '---- g_rhs-aft-get_prelim_fock -------- START' 342 call ga_print(g_rhs) 343 call ga_print(g_rhs0) 344 call ga_print(g_d1) 345 if (ga_nodeid().eq.0) 346 & write(*,*) '---- g_rhs-aft-get_prelim_fock -------- START' 347 endif 348c 349c Build "first order fock matrix" 350c 351 if (use_theory.eq.'dft') then 352 if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .true.)) 353 $ call errquit('hnd_giaox_zora: rtdb_put of xc_active failed', 354 & 0,RTDB_ERR) 355 if(.not. rtdb_put(rtdb,'fock_xc:calc_type', MT_INT, 1, 2)) 356 $ call errquit('hnd_giaox_zora: rtdb_put of calc_type failed', 357 & 0,RTDB_ERR) 358 if(.not. rtdb_put(rtdb,'fock_j:derfit', MT_LOG, 1, .false.)) 359 $ call errquit('hnd_giaox_zora: rtdb_put of j_derfit failed',0, 360 & RTDB_ERR) 361 endif 362 clo(1) = 3*npol*2 363 clo(2) = nbf 364 clo(3) = nbf 365 chi(1) = 1 366 chi(2) = -1 367 chi(3) = -1 368 if (.not.nga_create(MT_DBL,3,clo,'Fock matrix',chi,g_fock)) call 369 & errquit('hnd_giaox_zora: nga_create failed g_fock',0,GA_ERR) 370 call ga_zero(g_fock) 371c 372c Note: Just the exchange: jfac = 0.d0 (see above) 373c 374 if (.not.cam_exch) then 375 call shell_fock_build(geom, basis,0,3*npol*2, 376 $ jfac,kfac,tol2e, 377 & g_d1, ! input 378 & g_fock,! output 379 & .false.) 380 else 381 call shell_fock_build_cam(geom, basis,0,3*npol*2, 382 $ jfac,kfac,tol2e, 383 & g_d1, ! input 384 & g_fock,! output 385 & .false.) 386 end if 387c 388 if(use_theory.eq.'dft') then 389 if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, 0)) 390 $ call errquit('hnd_giaox_zora: rtdb_put failed',0,RTDB_ERR) 391 endif 392 if (npol.eq.1) then 393 nocc(1)=nclosed(1) 394 nocc(2)=0 395 else if (npol.eq.2) then 396 nocc(1)=nopen(1) 397 nocc(2)=nopen(2) 398 endif 399 if (debug_giaox.eq.1) then 400 if (ga_nodeid().eq.0) 401 & write(*,*) '---- g_rhs-bef-add_fock -------- START' 402 call ga_print(g_rhs) 403 call ga_print(g_fock) 404 if (ga_nodeid().eq.0) 405 & write(*,*) '---- g_rhs-bef-add_fock -------- END' 406 endif 407 call add_fock(g_rhs, ! out: g_rhs=g_rhs+g_fock 408 & g_fock,vectors, 409 & nbf,nmo,npol,nocc,nvirt) 410 411 if (debug_giaox.eq.1) then 412 if (ga_nodeid().eq.0) then 413 write(*,*) 'Aft add_fock()' 414 write(*,*) '---- g_rhs-aft-add_fock -------- START' 415 endif 416 call ga_print(g_rhs) 417 if (ga_nodeid().eq.0) 418 & write(*,*) '---- g_rhs-aft-add_fock -------- END' 419 endif 420 421c Cleanup of g_d1 and g_fock, not needed for now 422 if (.not.ga_destroy(g_d1)) call 423 & errquit('hnd_giaox_zora: ga_destroy failed g_d1',0,GA_ERR) 424 if (.not.ga_destroy(g_fock)) call 425 & errquit('hnd_giaox_zora: ga_destroy failed g_fock',0,GA_ERR) 426 427 if (debug_giaox.eq.1) then 428 write(*,*) '---- g_rhs-bef-add_H10 -------- START' 429 call ga_print(g_rhs) 430 write(*,*) '---- g_rhs-bef-add_H10 -------- END' 431 endif 432c ------ WARNING: skip new_giao_2d() ---- START 433c if (ga_nodeid().eq.0) 434c & write(*,*) 'WARNING: SKIP add_H10() ...' 435c goto 178 436c ------ WARNING: skip new_giao_2d() ---- END 437 438 if (switch_nmrcs_analysis) then 439 call add_H10_debug( 440 & g_rhs, ! out: accumulated rhs expression 441 & g_rhs_1e, ! out : Fji^{k,1e} 442 & ga_Fji, ! in: Fock 1st-deriv without V (pot.) contrib. 443 & vectors, ! in: MO coeffs 444 & dbl_mb(k_xyz), ! in: Nuclear positions (x,y,z) 445 & nat_slc, ! in: nr. selected nuclei (atoms) 446 & basis, ! in: basis handle 447 & nbf, ! in: nr. basis functions 448 & nmo, ! in: nr. MOs (occ+virt) 449 & npol, ! in: nr. of polarizations 450 & nocc, ! in: nr. occ MOs 451 & nvirt, ! in: nr. virtual MOs 452 & do_zora, ! in: switch = .true. zora on 453 & rtdb) ! in: rtdb handle 454 else 455 call add_H10(g_rhs, ! out: ga_rhs(a,i)=ga_rhs(a,i)+H10(a,i) 456 & ga_Fji,vectors, 457 & dbl_mb(k_xyz),nat_slc,basis, 458 & nbf,nmo,npol,nocc,nvirt,do_zora,rtdb) 459 endif 460c178 continue 461 462 if (debug_giaox.eq.1) then 463 if (ga_nodeid().eq.0) then 464 write(*,*) 'Aft. add_H10()' 465 write(*,*) '---- g_rhs-aft-add_H10 -------- START' 466 endif 467 call ga_print(g_rhs) 468 if (ga_nodeid().eq.0) 469 & write(*,*) '---- g_rhs-aft-add_H10 -------- END' 470 endif 471c ------------ added-4 ------------ END 472c Remaining term is Perturbed (GIAO) two-electron term times Unperturbed density 473c Calculate Sum(r,s) D0(r,s) * G10(m,n,r,s) in AO basis 474 alo(1) = -1 475 alo(2) = -1 476 alo(3) = 1 477 ahi(1) = nbf 478 ahi(2) = nbf 479 ahi(3) = 3*npol 480 if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix',alo,g_fock)) call 481 & errquit('hnd_giaox_zora: nga_create failed g_fock',0,GA_ERR) 482 call ga_zero(g_fock) 483 if (switch_nmrcs_analysis) then 484 if (ga_nodeid().eq.0) 485 & write(*,*) 'ENTER-define g_fock_Coul,g_fock_Exch' 486c -- TEST: create g_fock_Coul,g_fock_Exch ----------- START 487 if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix', 488 & alo,g_fock_Coul)) call 489 & errquit('hnd_giaox_zora: nga_create failed g_fockJ',0,GA_ERR) 490 call ga_zero(g_fock_Coul) 491 if (.not.nga_create(MT_DBL,3,ahi,'Fock matrix', 492 & alo,g_fock_Exch)) call 493 & errquit('hnd_giaox_zora: nga_create failed g_fockX',0,GA_ERR) 494 call ga_zero(g_fock_Exch) 495c -- TEST: create g_fock_Coul,g_fock_Exch ----------- END 496 endif 497 if(use_theory.eq.'dft') then 498 ifld = 4 499 if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, ifld)) 500 $ call errquit('hnd_giaox_zora: rtdb_put failed',0,RTDB_ERR) 501 endif 502c +++++++++++++++++++++++++++++++++++++ 503c ++++++++ calling new_giao_2e() ++++++ 504 if (debug_giaox.eq.1) then 505 if (ga_nodeid().eq.0) 506 & write(*,*) '---- g_dens bef new_giao_2e -------- START' 507 do cdens=1,npol 508 call ga_print(g_dens(cdens)) 509 enddo 510 if (ga_nodeid().eq.0) 511 & write(*,*) '---- g_dens bef new_giao_2e -------- END' 512 endif 513 call ga_sync() 514c ------ WARNING: skip new_giao_2d() ---- START 515c if (ga_nodeid().eq.0) 516c & write(*,*) 'WARNING: SKIP new_giao_2e() ...' 517c goto 176 518c ------ WARNING: skip new_giao_2d() ---- END 519 520 if (switch_nmrcs_analysis) then 521c +++++ Test: getting Coulomb and Exchange contrib separate --- START 522 call new_giao_2e_JK(geom,basis,nbf,tol2e, 523 & g_dens, ! in: e-denstiy 524 & g_fock, ! out: fock matrix 525 & g_fock_Coul, 526 & g_fock_Exch, 527 & xfac, 528 & npol) 529c +++++ Test: getting Coulomb and Exchange contrib separate --- END 530 else ! default 531 call new_giao_2e(geom,basis,nbf,tol2e, 532 & g_dens, ! in: e-denstiy 533 & g_fock, ! out: fock matrix 534 & xfac, 535 & npol) 536 endif 537 176 continue 538 539c ++++++++ calling new_giao_2e() ++++++ 540c +++++++++++++++++++++++++++++++++++++ 541 542 if (debug_giaox.eq.1) then 543 if (ga_nodeid().eq.0) 544 & write(*,*) '---- g_fock-aft-new_giao -------- START' 545 call ga_print(g_fock) 546 if (ga_nodeid().eq.0) 547 & write(*,*) '---- g_fock-aft-new_giao -------- END' 548 endif 549 550 if(use_theory.eq.'dft') then 551 ifld = 0 552 if (.not. rtdb_put(rtdb, 'fock_xc:calc_type', mt_int, 1, ifld)) 553 $ call errquit('hnd_giaox_zora: rtdb_put failed',0,RTDB_ERR) 554 if(.not. rtdb_put(rtdb,'bgj:xc_active', MT_LOG, 1, .false.)) 555 $ call errquit('hnd_giaox_zora: rtdb_put of xc_active failed', 556 & 0,RTDB_ERR) 557 endif 558 if (debug_giaox.eq.1) then 559 write(*,114) npol,nocc(1),nocc(2), 560 & nclosed(1),nclosed(2), 561 & nvirt(1),nvirt(2),scftyp 562 114 format('BEF. giao_aotomo: npol=',i3, 563 & ' nocc=(',i3,',',i3,') ', 564 & 'nclos=(',i3,',',i3,') ', 565 & 'nvirt=(',i3,',',i3,') scftyp=',a,')') 566 endif 567c Transform to MO basis and add to right-hand-side 568 call giao_aotomo(g_fock,vectors,nocc,nvirt,npol,3,nbf) 569 if (switch_nmrcs_analysis) then 570c ----- TEST: AO2MO for g_fock_Coul,g_fock_Exch ------ START 571 call giao_aotomo(g_fock_Coul,vectors,nocc,nvirt,npol,3,nbf) 572 call giao_aotomo(g_fock_Exch,vectors,nocc,nvirt,npol,3,nbf) 573c ----- TEST: AO2MO for g_fock_Coul,g_fock_Exch ------ END 574 endif 575 if (debug_giaox.eq.1) then 576 write(*,*) 'Aft add_fock1-giao-aotomo()' 577 write(*,*) 'BEF. giao_aotomo: npol=' 578 write(*,*) '---- g_fock-aft-giao_aotomo -------- START' 579 call ga_print(g_fock) 580 write(*,*) '---- g_fock-aft-giao_aotomo -------- END' 581 endif 582 583 if (switch_nmrcs_analysis) then 584c ------- Store g_rhs without Coulomb and Exchange ---- START 585 call ga_copy(g_rhs,g_rhs_noJK) 586 do ispin=1,npol 587c ------ definitions for g_rhs -------- START 588 disp = nocc(1)*nvirt(1)*(ispin-1) 589 shift=3*(ispin-1) 590 blo(1) = disp+1 591 bhi(1) = disp+nocc(ispin)*nvirt(ispin) 592 blo(2) = 1 593 bhi(2) = 3 594c ------ definitions for g_rhs -------- END 595 alo(1) = nocc(ispin)+1 596 ahi(1) = nmo 597 alo(2) = 1 598 ahi(2) = nocc(ispin) 599 alo(3) = shift+1 600 ahi(3) = shift+3 601 if (npol.eq.1) then 602 call nga_scale_patch(g_rhs_noJK,blo,bhi,-4.0d0) 603 call nga_scale_patch(g_rhs_1e ,blo,bhi,-4.0d0) 604 call nga_scale_patch(g_rhs_eSji,blo,bhi,-4.0d0) 605 else if (npol.eq.2) then 606 call nga_scale_patch(g_rhs_noJK,blo,bhi,-2.0d0) 607 call nga_scale_patch(g_rhs_1e ,blo,bhi,-2.0d0) 608 call nga_scale_patch(g_rhs_eSji,blo,bhi,-2.0d0) 609 endif 610 enddo ! end-loop-ispin 611c ------- Store g_rhs without Coulomb and Exchange ---- END 612 endif 613 614 call add_fock1(g_rhs, ! out: accumulated rhs expression 615 & g_fock,! in 616 & nmo,npol,nocc,nvirt) 617 618 if (switch_nmrcs_analysis) then 619 call add_fock1(g_rhs_Coul, ! out: accumulated rhs expression 620 & g_fock_Coul,! in 621 & nmo,npol,nocc,nvirt) 622 call add_fock1(g_rhs_Exch, ! out: accumulated rhs expression 623 & g_fock_Exch,! in 624 & nmo,npol,nocc,nvirt) 625 if (.not. ga_destroy(g_fock_Coul)) call errquit( 626 & 'hnd_giaox_zora: ga_destroy failed g_fock_Coul',0, GA_ERR) 627 if (.not. ga_destroy(g_fock_Exch)) call errquit( 628 & 'hnd_giaox_zora: ga_destroy failed g_fock_Exch',0, GA_ERR) 629 endif 630 if (debug_giaox.eq.1) then 631 write(*,*) 'Aft add_fock1()' 632 write(*,*) '---- g_rhs-aft-add_fock1 -------- START' 633 call ga_print(g_rhs) 634 write(*,*) '---- g_rhs-aft-add_fock1 -------- END' 635 endif 636 if (.not.ga_destroy(g_fock)) call 637 & errquit('hnd_giaox_zora: ga_destroy failed g_fock',0,GA_ERR) 638 639 if (switch_skip_cphf) then 640 if (switch_nmrcs_analysis) then 641 if (ga_nodeid().eq.0) 642 & write(*,*) 'WARNING : Using skip_cphf_JK ...' 643 call skip_cphf_JK( 644 & g_rhs, ! IN/OUT 645 & g_rhs_Coul, 646 & g_rhs_Exch, 647 & g_rhs_noJK, 648 & g_rhs_1e, 649 & g_rhs_eSji, 650 & dbl_mb(k_eval), ! IN: energy eigenvalues 651 & nocc, 652 & nvirt, 653 & nbf, ! FA-04-28-12: replacing nmo by nbf 654 & npol) 655 else 656 if (ga_nodeid().eq.0) 657 & write(*,*) 'WARNING : Using skip_cphf ...' 658 call skip_cphf(g_rhs, ! IN/OUT 659 & dbl_mb(k_eval), ! IN: energy eigenvalues 660 & nocc, 661 & nvirt, 662 & nbf, ! FA-04-28-12: replacing nmo by nbf 663 & npol) 664 endif 665 endif 666 667c -------free allocated memory -------------- START 668 if (.not.ma_pop_stack(l_eval)) call 669 & errquit('hnd_giaox_zora: ma_pop_stack failed k_eval',0,MA_ERR) 670 if (.not.ma_pop_stack(l_occ)) call 671 & errquit('hnd_giaox_zora: ma_pop_stack failed k_occ',0,MA_ERR) 672c -------free allocated memory -------------- END 673 if (debug_giaox.eq.1) then 674 if (ga_nodeid().eq.0) 675 & write(*,*) '---- Reading INPUT-cphf: g_rhs -------- START' 676 call ga_print(g_rhs) 677 if (ga_nodeid().eq.0) 678 & write(*,*) '---- Reading INPUT-cphf: g_rhs -------- END' 679 endif 680 if(.not.rtdb_get(rtdb,'zora:skip_cphf_ev_shield', 681 & mt_log,1,skip_cphf_ev_shield)) 682 & skip_cphf_ev_shield = .false. 683 if (.not.(switch_skip_cphf) .and. 684 & .not.(skip_cphf_ev_shield)) then 685c if (ga_nodeid().eq.0) 686c & write(*,*) 'COMPUTE cphf shield data ...' 687c Write ga_rhs to disk 688 call cphf_fname('cphf_rhs',cphf_rhs) 689 call cphf_fname('cphf_sol',cphf_sol) 690 if(.not.file_write_ga(cphf_rhs,g_rhs)) call errquit 691 $ ('hnd_giaox_zora: could not write cphf_rhs',0, DISK_ERR) 692 call schwarz_tidy() 693 call int_terminate() 694c 695c Call the CPHF routine 696c 697c We do need to tell the CPHF that the density is skew symmetric. 698c Done via rtdb, put cphf:skew .false. on rtdb and later remove it. 699 700 if (.not. rtdb_put(rtdb, 'cphf:skew', mt_log, 1,.false.)) call 701 $ errquit('hnd_giaox_zora: failed to write skew ', 0, RTDB_ERR) 702 if (.not.cphf2(rtdb)) call errquit 703 $ ('hnd_giaox_zora: failure in cphf ',0, RTDB_ERR) 704 if (.not. rtdb_delete(rtdb, 'cphf:skew')) call 705 $ errquit('hnd_giaox_zora: rtdb_delete failed ', 0, RTDB_ERR) 706c 707c Occ-virt blocks are the solution pieces of the CPHF 708c Read solution vector from disk and put solutions in U matrices 709 call ga_zero(g_rhs) 710 if(.not.file_read_ga(cphf_sol,g_rhs)) call errquit 711 $ ('hnd_giaox_zora: could not read cphf_rhs',0, DISK_ERR) 712c ----- write CPHF data to file ----------------- START 713c Note.- lbl_cphfhyp defined in zora.fh 714 call util_file_name(lbl_cphfshield, 715 & .false.,.false.,zorafilename) 716c if (ga_nodeid().eq.0) 717c & write(*,*) '---------- g_rhs0-cs --------- START' 718c call ga_print(g_rhs0) 719c if (ga_nodeid().eq.0) 720c & write(*,*) '---------- g_rhs0-cs --------- END' 721c if (ga_nodeid().eq.0) 722c & write(*,*) '---------- g_rhs-cs --------- START' 723c call ga_print(g_rhs) 724c if (ga_nodeid().eq.0) 725c & write(*,*) '---------- g_rhs-cs --------- END' 726 call dft_zoraCPHF_write( 727 & zorafilename, ! in: filename 728 & npol, ! in: nr polarization 729 & nocc, ! in: nr occupied MOs 730 & nvirt, ! in: nr virtual MOs 731 & nbf, ! in: nr basis functions 732 & vectors, ! in: MOs 733 & g_rhs0, ! in: (ntot,3) GA matrix 734 & g_rhs) ! in: (nocc*nvirt,3) GA matrix 735c ----- write CPHF data to file ----------------- END 736 else 737 if (ga_nodeid().eq.0) 738 & write(*,*) 'WARNING: SKIP cphf ...' 739 if (skip_cphf_ev_shield) then 740 call ga_zero(g_rhs) 741 do i=1,npol 742 call ga_zero(vectors(i)) 743 enddo 744 if (ga_nodeid().eq.0) 745 & write(*,*) 'READ cphf shieldings data from file ...' 746 call util_file_name(lbl_cphfshield,.false.,.false.,zorafilename) 747 call dft_zoraCPHF_read( 748 & zorafilename, ! in: filename 749 & npol, ! in: nr polarization 750 & nocc, ! in: nr occupied MOs 751 & nvirt, ! in: nr virtual MOs 752 & nbf, ! in: nr basis functions 753 & vectors, ! out: MOs 754 & g_rhs0, ! out: (ntot,3) GA matrix 755 & g_rhs) ! out: (nocc*nvirt,3) GA matrix 756c if (ga_nodeid().eq.0) 757c & write(*,*) '---------- g_rhs-cs-read --------- START' 758c call ga_print(g_rhs) 759c if (ga_nodeid().eq.0) 760c & write(*,*) '---------- g_rhs-cs-read --------- END' 761 endif 762 endif ! end-if-switch-skip_cphf 763 764 if (debug_giaox.eq.1) then 765 if (ga_nodeid().eq.0) 766 & write(*,*) '---- Reading sol-cphf: g_rhs -------- START' 767 call ga_print(g_rhs) 768 if (ga_nodeid().eq.0) 769 & write(*,*) '---- Reading sol-cphf: g_rhs -------- END' 770 endif 771 type_NMR=1 ! =1,2,3=shieldings,hyperfine,gshift 772 if (switch_nmrcs_analysis) then 773 call get_P10_JK( 774 & g_d1_oo, ! out: Perturbed density matrix occ-occ contrib 775 & g_d1_ov, ! out: Perturbed density matrix occ-virt contrib 776 & type_NMR, 777 & g_d1_ov_Coul,! out: Perturbed (spin)density matrix with Coulomb contrib 778 & g_d1_ov_Exch,! out: Perturbed (spin)density matrix with Exchange contrib (fock_xc) 779 & g_d1_ov_noJK,! out: Perturbed (spin)density matrix remaining terms not J or K 780 & g_d1_ov_1e, ! out: Perturbed (spin)density matrix 1e-operator contrib 781 & g_d1_ov_eSji,! out: Perturbed (spin)density matrix perturbed overlap matrix 782 & g_rhs, ! in: accumulated rhs expression 783 & g_rhs_Coul, 784 & g_rhs_Exch, 785 & g_rhs_noJK, 786 & g_rhs_1e, 787 & g_rhs_eSji, 788 & g_rhs0, ! in: from get_prelim_fock() 789 & vectors,g_CiFull, 790 & nbf,nmo,npol,nocc,nvirt, 791 & do_zora,do_NonRel,not_zora_scale, 792 & lbl_nlmogshift, ! in: for g-shift NLMO analysis 793 & lbl_nlmoshield, ! in: for shield NLMO analysis 794 & rtdb) 795 if (.not. ga_destroy(g_rhs_Coul)) call errquit( 796 & 'hnd_giaox_zora: ga_destroy failed g_rhs_Exch',0, GA_ERR) 797 if (.not. ga_destroy(g_rhs_Exch)) call errquit( 798 & 'hnd_giaox_zora: ga_destroy failed g_rhs_Exch',0, GA_ERR) 799 if (.not. ga_destroy(g_rhs_noJK)) call errquit( 800 & 'hnd_giaox_zora: ga_destroy failed g_rhs_noJK',0, GA_ERR) 801 if (.not. ga_destroy(g_rhs_1e)) call errquit( 802 & 'hnd_giaox_zora: ga_destroy failed g_rhs_1e',0, GA_ERR) 803 if (.not. ga_destroy(g_rhs_eSji)) call errquit( 804 & 'hnd_giaox_zora: ga_destroy failed g_rhs_1e',0, GA_ERR) 805 else ! default 806c call get_P10( g_d1, ! out: Perturbed density matrix 807c & g_rhs, ! in: accumulated rhs expression 808c & g_rhs0, ! in: from get_prelim_fock() 809c & vectors,g_CiFull, 810c & nbf,nmo,npol,nocc,nvirt,rtdb) 811 call get_P10_1( 812 & g_d1_oo, ! out: Perturbed density matrix occ-occ contrib 813 & g_d1_ov, ! out: Perturbed density matrix occ-virt contrib 814 & type_NMR, ! in : =1,2,3=shieldings,hyperfine,gshift 815 & g_rhs, ! in: accumulated rhs expression 816 & g_rhs0, ! in: from get_prelim_fock() 817 & vectors,g_CiFull, 818 & nbf,nmo,npol,nocc,nvirt, 819 & do_zora,do_NonRel,not_zora_scale, 820 & lbl_nlmogshift, ! in : label for g-shift NLMO analysis 821 & lbl_nlmoshield, ! in : label for shield NLMO analysis 822 & rtdb) 823 endif 824 825 if (.not.ga_destroy(g_rhs)) call 826 & errquit('hnd_giaox_zora: ga_destroy failed g_rhs',0,GA_ERR) 827 if (.not.ga_destroy(g_rhs0)) call 828 & errquit('hnd_giaox_zora: ga_destroy failed g_rhs0',0,GA_ERR) 829 830c Now we have in g_d1(nmo,nmo,3) the derivative densities and 831c hence we can calculate the contributions to the shielding tensor 832 if (.not. ma_push_get(mt_dbl,9*nat_slc,'sh para',l_para,k_para)) 833 & call errquit('hnd_giaox_zora: ma_push_get failed k_para', 834 & 0,MA_ERR) 835 if (.not. ma_push_get(mt_dbl,9*nat_slc,'sh dia',l_dia,k_dia)) call 836 & errquit('hnd_giaox_zora: ma_push_get failed k_dia',0,MA_ERR) 837c Before we start getting the integrals we need to reinitialize the 838c integrals. They were terminated by the cphf. 839 call int_init(rtdb,1,basis) 840 call schwarz_init(geom,basis) 841 call hnd_giao_init(basis,1) 842 843c ------- TEST: get_par() get_dia()------ START 844 if (switch_nmrcs_analysis) then 845 npar_analysis=9 846 else 847 npar_analysis=4 848 endif 849 if (.not. ma_push_get(mt_dbl, 850 & npar_analysis*3*3*nat_slc,'sh pararr', 851 & l_pararr,k_pararr)) 852 & call errquit('hnd_giaox_zora: ma_push_get failed k_pararr', 853 & 0,MA_ERR) 854 if (.not. ma_push_get(mt_dbl, 855 & 3*3*nat_slc,'sh fukui', 856 & l_fukui,k_fukui)) 857 & call errquit('hnd_giaox_zora: ma_push_get failed l_fukui', 858 & 0,MA_ERR) 859 call ga_sync() 860 if (switch_nmrcs_analysis) then 861 call get_par_JK( 862 & dbl_mb(k_para), ! OUT : paramagnetic tensor 863 & ga_para1, ! IN : paramagnetic tensor (GA) 864 & ga_h01_num, ! IN : H01 matrix 865 & dbl_mb(k_pararr),! OUT : par tensor gauge,OO,OV 866 & g_d1_oo, ! IN : Perturbed density matrix (OO part) 867 & g_d1_ov, ! IN : Perturbed density matrix (OV part) 868 & g_d1_ov_Coul, ! IN : Perturbed density matrix (OV part-Coul only) 869 & g_d1_ov_Exch, ! IN : Perturbed density matrix (OV part-Exch only) 870 & g_d1_ov_noJK, ! IN : Perturbed density matrix (OV part-all -(Coul,Exch)) 871 & g_d1_ov_1e, ! IN : Perturbed density matrix (OV part-1e contrib) 872 & g_d1_ov_eSji, ! IN : Perturbed density matrix (OV part-eSji contrib) 873 & g_dens, 874 & npol, ! IN : nr. polarizations 875 & dbl_mb(k_fukui), 876 & basis, ! IN : basis handle 877 & nbf, ! IN : nr. basis functions 878 & dbl_mb(k_xyz), ! IN : Nuclear positions (x,y,z) 879 & nat_slc, ! IN : nr. atoms selected 880 & nat, ! IN : nr. atoms total 881 & geom, 882 & rtdb, 883 & oskel) 884 if (.not. ga_destroy(g_d1_ov_Coul)) call errquit( 885 & 'hnd_giaox_zora: ga_destroy failed g_d1_ov_Coul',0, GA_ERR) 886 if (.not. ga_destroy(g_d1_ov_Exch)) call errquit( 887 & 'hnd_giaox_zora: ga_destroy failed g_d1_ov_Exch',0, GA_ERR) 888 if (.not. ga_destroy(g_d1_ov_noJK)) call errquit( 889 & 'hnd_giaox_zora: ga_destroy failed g_d1_ov_noJK',0, GA_ERR) 890 if (.not. ga_destroy(g_d1_ov_1e)) call errquit( 891 & 'hnd_giaox_zora: ga_destroy failed g_d1_ov_noJK',0, GA_ERR) 892 if (.not. ga_destroy(g_d1_ov_eSji)) call errquit( 893 & 'hnd_giaox_zora: ga_destroy failed g_d1_ov_noJK',0, GA_ERR) 894 else 895 call get_par( 896 & dbl_mb(k_para), ! OUT : paramagnetic tensor 897 & ga_para1, ! IN : paramagnetic tensor (GA) 898 & ga_h01_num, ! IN : H01 matrix 899 & dbl_mb(k_pararr),! OUT : par tensor gauge,OO,OV 900 & g_d1_oo, ! IN : Perturbed density matrix (OO part) 901 & g_d1_ov, ! IN : Perturbed density matrix (OV part) 902 & g_dens, ! IN : e-density 903 & npol, ! IN : nr. polarizations 904 & dbl_mb(k_fukui), ! OUT : Fukui term 905 & basis, ! IN : basis handle 906 & nbf, ! IN : nr. basis functions 907 & dbl_mb(k_xyz), ! IN : Nuclear positions (x,y,z) 908 & nat_slc, ! IN : nr. atoms selected 909 & nat, ! IN : nr. atoms total 910 & geom, 911 & rtdb, 912 & oskel) ! IN : = .false. 913 endif 914 915 call get_dia(dbl_mb(k_dia), ! OUT: dia tensor 916 & ga_dia, ! IN : ga dia tensor 917 & basis, ! IN : basis handle 918 & g_dens, ! IN : e-density 919 & nbf, ! IN : nr. basis functions 920 & npol, ! IN : nr. polarizations 921 & dbl_mb(k_fukui), ! IN : Fukui term 922 & dbl_mb(k_xyz), ! IN : nuclear positions 923 & nat_slc, ! IN : selected atoms 924 & rtdb, ! IN : rtdb handle 925 & oskel) ! IN : = .false 926c ------- TEST: get_par() get_dia()------ END 927c +++++++++ print-total-pardia-transferred +++ START 928 debug_cs=.false. 929 if (debug_cs) then 930 if (ga_nodeid().eq.0) then 931 ic=1 932 do iatom = 1, nat_slc 933 do ix = 1, 3 934 do iy = 1, 3 935 cbuf=k_pararr+ 936 & 3*3*npar_analysis*(iatom-1)+ 937 & 3*npar_analysis*(iy-1)+npar_analysis*(ix-1)-1 938 if (switch_nmrcs_analysis) then !-- START-if-switch_gshift_analysis 939 write(*,179) ix,iy,iatom, 940 & dbl_mb(k_dia+ic-1), 941 & dbl_mb(cbuf+1),dbl_mb(cbuf+2), 942 & dbl_mb(cbuf+3), 943 & dbl_mb(cbuf+5),dbl_mb(cbuf+6), 944 & dbl_mb(cbuf+7), 945 & dbl_mb(cbuf+8),dbl_mb(cbuf+9), 946 & dbl_mb(cbuf+4) 947 179 format('NW:(dia,gauge,OO,OV,', 948 & 'OV_Coul,OV_Exch,OV_nJK,', 949 & 'OV_1e,OV_eSji,' 950 & 'Totpar)(',i1,',',i1,',',i1,')=(', 951 & f18.6,' ',f18.6,' ',f18.6,' ', 952 & f18.6,' ',f18.6,' ', 953 & f18.6,' ',f18.6,' ', 954 & f18.6,' ',f18.6,' ', 955 & f18.6,' )') 956 else 957 write(*,19) ix,iy,iatom, 958 & dbl_mb(k_dia+ic-1), 959 & dbl_mb(cbuf+1),dbl_mb(cbuf+2), 960 & dbl_mb(cbuf+3),dbl_mb(cbuf+4), 961 & dbl_mb(k_dia+ic-1)+dbl_mb(cbuf+4) 962 19 format('NW:(dia,gauge,OO,OV,Totpar,dia+par)(', 963 & i1,',',i1,',',i1,')=(', 964 & f18.6,' ',f18.6,' ',f18.6,' ',f18.6,' ', 965 & f18.6,' ',f18.6,' )') 966 endif 967 ic=ic+1 968 enddo ! end-loop-iy 969 enddo ! end-loop-ix 970 enddo ! end-loop-iatom 971 endif ! end-if-ga_nodeid-eq-0 972 endif ! end-if-debug_cs 973c +++++++++ print-total-dia-transferred +++ END 974 if (.not.ma_pop_stack(l_fukui)) call 975 & errquit('hnd_giaox_zora: ma_pop_stack failed k_fukui', 976 & 0,MA_ERR) 977 if (.not.ma_pop_stack(l_pararr)) call 978 & errquit('hnd_giaox_zora: ma_pop_stack failed k_para',0,MA_ERR) 979 do ispin=1,ndens 980 if (.not.ga_destroy(g_dens(ispin))) call 981 & errquit('hnd_giaox_zora: ga_destroy failed g_dens', 982 & 0,GA_ERR) 983 enddo 984 985c --------- allocating array for NLMO analysis --------------- START 986 shldfile=0 ! not doing NLMO analysis by default 987 status=rtdb_get(rtdb,'prop:shldfile',mt_int,1,shldfile) ! for NLMO analysis 988 if (shldfile.eq.1) then 989 if (.not. ma_alloc_get(mt_dbl,nat_slc*9,'tvec',l_tvec,k_tvec)) 990 & call 991 & errquit('hnd_giaox_zora: ma_push_get failed tvec',0,MA_ERR) 992 endif ! end-if-shldfile 993c --------- allocating array for NLMO analysis --------------- END 994c 995c Print out tensor information, and write to Ecce file if necessary 996 status = rtdb_parallel(.false.) 997 if (ga_nodeid().gt.0) goto 300 998 acc_vec=0 ! For NMLO analysis 999 call ecce_print_module_entry('nmr') 1000 do iatom = 1, nat_slc 1001 ioff = (iatom-1)*9 1002 if (.not. geom_cent_get(geom, int_mb(k_AtNr+iatom-1), tag, 1003 & dbl_mb(k_xyz), dbl_mb(k_zan))) call 1004 & errquit('hnd_giaox_zora: geom_cent_tag failed',0,GEOM_ERR) 1005 if (.not. geom_tag_to_element(tag, symbol, element, atn)) then 1006 if (.not. inp_compare(0,tag(1:2),'bq')) call ! check for bq as a fall back 1007 & errquit('hnd_giaox_zora: geom_tag_to_element failed', 1008 & 0,GEOM_ERR) 1009 endif 1010c 1011c Print tensor pieces and sum for total shielding tensor 1012 if (ga_nodeid().eq.0) then 1013 write(luout,9700) iatom,symbol ! Showing original atm nr. 1014 write(luout,9800) (dbl_mb(k_dia+ioff+ix-1),ix=1,9) 1015 write(luout,9801) (dbl_mb(k_para+ioff+ix-1),ix=1,9) 1016 endif 1017 do ix = 0, 8 1018 dbl_mb(k_para+ioff+ix) = dbl_mb(k_dia +ioff+ix) + 1019 & dbl_mb(k_para+ioff+ix) 1020 enddo 1021c 1022c Print total shielding tensor 1023c 1024 if (ga_nodeid().eq.0) then 1025 write(luout,9802) (dbl_mb(k_para+ioff+ix-1),ix=1,9) 1026c 1027c Diagonalize total tensor 1028c Order in a: xx xy yy xz yz zz 1029 a(1) = dbl_mb(k_para+ioff) 1030 a(2) = dbl_mb(k_para+ioff+1) 1031 a(3) = dbl_mb(k_para+ioff+4) 1032 a(4) = dbl_mb(k_para+ioff+2) 1033 a(5) = dbl_mb(k_para+ioff+5) 1034 a(6) = dbl_mb(k_para+ioff+8) 1035 ij = 0 1036 do 241 i = 1, 3 1037 do 241 j = 1, i 1038 ij = ij + 1 1039 axs(i,j) = a(ij) 1040 axs(j,i) = a(ij) 1041 241 continue 1042 call hnd_diag(axs,eig,3,.true.,.true.) 1043 isotr =(eig(1) + eig(2) + eig(3))/3.0d0 1044 aniso = eig(1) -(eig(2) + eig(3))/2.0d0 1045c ++++++++++ get eigenvectors for NMLO analysis +++++ START 1046 shldfile=0 ! not doing NLMO analysis by default 1047 status=rtdb_get(rtdb,'prop:shldfile',mt_int,1,shldfile) ! for NLMO analysis 1048 if (shldfile.eq.1) then 1049 do i1=1,3 1050 do j1=1,3 1051 dbl_mb(k_tvec+acc_vec)=axs(i1,j1) 1052 acc_vec=acc_vec+1 1053 enddo 1054 enddo 1055 endif 1056c ++++++++++ get eigenvectors for NMLO analysis +++++ END 1057 write(luout,9987) isotr,aniso 1058 write(luout,9986) (ix,ix=1,3) 1059 write(luout,9985) (eig(ix),ix=1,3) 1060 do iy=1,3 1061 write(luout,9983) iy,(axs(iy,ix),ix=1,3) 1062 enddo 1063 write(luout,'(//)') 1064c 1065c Print Ecce information 1066c 1067 call ecce_print1_char('atom name',symbol,1) 1068 call ecce_print2('shielding tensor',MT_DBL, 1069 & dbl_mb(k_para+ioff),3,3,3) 1070 call ecce_print1('shielding isotropic' ,MT_DBL,isotr,1) 1071 call ecce_print1('shielding anisotropy' ,MT_DBL,aniso,1) 1072 call ecce_print1('shielding eigenvalues' ,MT_DBL,eig,3) 1073 call ecce_print2('shielding eigenvectors',MT_DBL,axs, 1074 & 3,3,3) 1075 endif 1076 enddo 1077 call ecce_print_module_exit('nmr','ok') 1078300 call ga_sync() 1079 1080 status = rtdb_parallel(.true.) 1081 shldfile=0 ! not doing NLMO analysis by default 1082 status=rtdb_get(rtdb,'prop:shldfile',mt_int,1,shldfile) ! for NLMO analysis 1083 if (shldfile.eq.1) then ! ------- hypfile-if++++ START 1084 if (.not. ga_create(mt_dbl,1,nat_slc*9, 1085 & 'munu4nbo: g_tvec',0,0,g_tvec)) 1086 $ call errquit('hnd_giaox_zora: g_tvec', 0,GA_ERR) 1087 call ga_dgop(msg_efgs_col,dbl_mb(k_tvec),nat_slc*9,'+') 1088 call ga_put(g_tvec,1,1,1,nat_slc*9,dbl_mb(k_tvec),1) 1089 call create_munu4nbo_shield( 1090 & rtdb, 1091 & g_tvec, 1092 & nat_slc,int_mb(k_AtNr), 1093 & basis,npol,nocc,nvirt,nmo) 1094 if (.not. ga_destroy(g_tvec)) call errquit( ! destroy GA 1095 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1096 if (.not.ma_free_heap(l_tvec)) call 1097 & errquit('hnd_giaox_zora: ma_free_heap l_tvec',0,MA_ERR) 1098 endif ! ------------------------ hypfile-if++++ END 1099c ---- Destroy stored ga arrays ------ START 1100 if (.not. ga_destroy(g_d1_oo)) call errquit( 1101 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1102 if (.not. ga_destroy(g_d1_ov)) call errquit( 1103 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1104 if (do_zora) then 1105 if (.not. ga_destroy(ga_dia)) call errquit( 1106 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1107 if (.not. ga_destroy(ga_para1)) call errquit( 1108 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1109 if (.not. ga_destroy(ga_h01_num)) call errquit( 1110 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1111 if (.not. ga_destroy(ga_Fji)) call errquit( 1112 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1113 if (.not. ga_destroy(g_AtNr1)) call errquit( 1114 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 1115 endif 1116c ---- Destroy stored ga arrays ------ END 1117c 1118c Clean up all remaining memory 1119 if (.not.ma_pop_stack(l_dia)) call 1120 & errquit('hnd_giaox_zora: ma_pop_stack failed k_dia',0,MA_ERR) 1121 if (.not.ma_pop_stack(l_para)) call 1122 & errquit('hnd_giaox_zora: ma_pop_stack failed k_para',0,MA_ERR) 1123 do i=1,npol 1124 911 if (.not.ga_destroy(vectors(i))) call 1125 & errquit('hnd_giaox_zora: ga_destroy failed vectors', 1126 & 0,GA_ERR) 1127 enddo 1128c if (.not.ga_destroy(vectors_scl(1))) call 1129c & errquit('giao_aotomo: ga_destroy failed vectors',0,GA_ERR) 1130 if (.not.ma_pop_stack(l_zan)) call 1131 & errquit('hnd_giaox_zora: ma_pop_stack failed k_zan',0,MA_ERR) 1132 if (.not.ma_pop_stack(l_xyz)) call 1133 & errquit('hnd_giaox_zora: ma_pop_stack failed k_xyz',0,MA_ERR) 1134 if (.not.ma_free_heap(l_AtNr)) call 1135 & errquit('hnd_giaox_zora: ma_free_heap l_AtNr',0,MA_ERR) 1136 call schwarz_tidy() 1137 call int_terminate() 1138 return 1139 7000 format(/,10x,'NMR shielding cannot be calculated for UHF', 1140 1 ' or ROHF wave functions: use RHF') 1141 9700 format(6x,'Atom: ',i4,2x,a2) 1142 9800 format(8x,'Diamagnetic',/,3(3F12.4,/)) 1143 9801 format(8x,'Paramagnetic',/,3(3F12.4,/)) 1144 9802 format(8x,'Total Shielding Tensor',/,3(3F12.4,/)) 1145 9983 format(6x,i1,3x,3f12.4) 1146 9985 format(10x,3f12.4,/) 1147 9986 format(10x,'Principal Components and Axis System',/,10x, 1148 1 3(7x,i1,4x)) 1149 9987 format(10x,' isotropic = ',f12.4,/, 1150 1 10x,'anisotropy = ',f12.4,/) 1151 9999 format( 1152 1 /,10x,41(1h-),/, 1153 2 10x,'Chemical Shielding Tensors (GIAO, in ppm)',/, 1154 3 10x,41(1h-),/) 1155 end 1156 1157 subroutine get_par(par, ! OUT : paramagnetic tensor 1158 & ga_para1, ! IN : paramagnetic tensor (GA) 1159 & ga_h01_num,! IN : H01 matrix 1160 & par_arr, ! OUT : par tensor gauge,OO,OV 1161 & g_d1_oo, ! IN : Perturbed density matrix (OO part) 1162 & g_d1_ov, ! IN : Perturbed density matrix (OV part) 1163 & g_dens, ! IN : e-density 1164 & npol, ! IN : nr. polarizations 1165 & Fukui, ! OUT : Fukui term 1166 & basis, ! IN : basis handle 1167 & nbf, ! IN : nr. basis functions 1168 & pos, ! IN : Nuclear positions (x,y,z) 1169 & natoms_slc,! IN : nr. atoms selected 1170 & natoms_tot,! IN : nr. atoms total 1171 & geom, 1172 & rtdb, 1173 & oskel) ! IN : =.false. 1174c Purpose : Assemble NMR Chemical Shielding: paramagnetic tensor 1175c Author : Fredy Aquino 1176c Date : 03-03-11, 11-23-12 (adding Fukui term) 1177c Note.- Adaptation of hnd_giaox.F to hold 1178c - Nonrel(HF or DFT) in unrestricted/restricted calculations. 1179c - Relativistic (zora) in unrestricted/restricted calculations. 1180c 1181c Date : 02-05-14: Fixed bug when selecting atoms 1182c BEFORE : call get_chi_centers_ga(g_R,basis,nbf,geom,natoms_slc) 1183c AFTER : call get_chi_centers_ga(g_R,basis,nbf,geom,natoms_tot) 1184 implicit none 1185c 1186#include "errquit.fh" 1187#include "global.fh" 1188#include "mafdecls.fh" 1189#include "msgids.fh" 1190#include "rtdb.fh" 1191#include "apiP.fh" 1192#include "prop.fh" 1193#include "bgj.fh" 1194#include "bas.fh" 1195#include "stdio.fh" 1196#include "zora.fh" 1197c 1198c Global arrays variables defined in zora.fh 1199c --> ga_dia,ga_para1,ga_Fji,ga_h01_num 1200 integer rtdb,geom 1201 integer basis 1202 integer nbf,natoms_slc,natoms_tot 1203 integer g_dens(3) ! IN : electronic density 1204 double precision Fukui(3,3,natoms_slc) ! used in analyt shield 1205 integer npol ! OUT: Used only in analytical shieldings 1206 integer alo(3), ahi(3), blo(3), bhi(3), 1207 & dlo(3), dhi(3) 1208 integer ld(2),cbuf,cbuf1,cbuf2 1209 integer i,dims(3),chunk(3), 1210 & g_R(3),ispin,a,b 1211 double precision pos(3*natoms_slc),val 1212 integer iatom,ix,iy,ixy,ind 1213 double precision val_oo,val_ov 1214 integer g_h01, ! local 1215 & ga_h01_num 1216 integer g_d1_oo,g_d1_ov ! input: Perturbed density matrics OO and OV 1217 ! OO, due to occupied-occupied MOs 1218 ! OV, due to occupied-viritual MOs 1219 integer ga_para1,g_t1,g_t2,g_t3 1220 integer debug_get_par 1221 double precision par(*) 1222 logical oskel 1223 double precision ppm,par_arr(*) 1224 integer ind_tmn(2,3) 1225 data ind_tmn / 2, 3, ! tmn=123 1226 & 3, 1, ! tmn=231 1227 & 1, 2 / ! tmn=312 1228 external get_chi_centers_ga 1229 data ppm /26.62566914d+00/ 1230 1231c Warning: If I activate debugging (=1) 1232c I could get trouble in final outputs 1233c the content is affected if I attempt 1234c to print out variables. (FA-11-24-11) 1235 1236 debug_get_par=0 1237 1238 do ixy = 1, 9*natoms_slc 1239 par(ixy) = 0.0d0 ! initialize 1240 enddo 1241 if (do_zora) then 1242 if (ga_nodeid().eq.0) 1243 & write(*,*) 'Calc. par tensor-> zora' 1244c +++++ Initialize dbl_mb(k_para) with ga_para1 ++++ START 1245c ---- STORE: g_para1 --> dbl_mb(k_para) 1246 alo(1)=1 1247 ahi(1)=3 1248 alo(2)=1 1249 ahi(2)=3 1250 alo(3)=1 1251 ahi(3)=natoms_slc 1252 ld(1)=3 1253 ld(2)=3 1254 call nga_get(ga_para1,alo,ahi,par(1),ld) 1255c +++++ Initialize dbl_mb(k_para) with g_para1 ++++ END 1256 alo(1) = 1 1257 ahi(1) = nbf 1258 alo(2) = 1 1259 ahi(2) = nbf 1260 blo(1) = 1 1261 bhi(1) = nbf 1262 blo(2) = 1 1263 bhi(2) = nbf 1264 ixy = 0 1265 blo(3) = 0 1266 bhi(3) = 0 1267 ind=0 1268 do iatom = 1, natoms_slc 1269 do iy = 1, 3 1270 blo(3) = blo(3) + 1 1271 bhi(3) = bhi(3) + 1 1272 do ix = 1, 3 1273 alo(3) = ix 1274 ahi(3) = ix 1275 ixy = ixy + 1 1276 val_oo = nga_ddot_patch(g_d1_oo ,'n',alo,ahi, 1277 & ga_h01_num,'n',blo,bhi) 1278 val_ov = nga_ddot_patch(g_d1_ov ,'n',alo,ahi, 1279 & ga_h01_num,'n',blo,bhi) 1280 1281 cbuf=9*(iatom-1)+3*(ix-1)+iy !-1 transpose 1282c ----- store in par_arr ------- START 1283 par_arr(ind+1)=par(cbuf) 1284 par_arr(ind+2)=val_oo*ppm 1285 par_arr(ind+3)=val_ov*ppm 1286 par_arr(ind+4)=par(cbuf)+ 1287 & val_oo*ppm+ 1288 & val_ov*ppm 1289 1290 if (debug_get_par.eq.1) then 1291 if (ga_nodeid().eq.0) then 1292 write(*,14) iatom,iy,ix, 1293 & alo(1),alo(2),alo(3), 1294 & ahi(1),ahi(2),ahi(3), 1295 & blo(1),blo(2),blo(3), 1296 & bhi(1),bhi(2),bhi(3), 1297 & par_arr(ind+1), 1298 & par_arr(ind+2),par_arr(ind+3), 1299 & par_arr(ind+4) 1300 14 format('(iatom,iy,ix)=(',i3,',',i3,',',i3,') ', 1301 & 'alo=(',i3,',',i3,',',i3,') ', 1302 & 'ahi=(',i3,',',i3,',',i3,') ', 1303 & 'blo=(',i3,',',i3,',',i3,') ', 1304 & 'bhi=(',i3,',',i3,',',i3,') ', 1305 & 'para=(',f15.8,',',f15.8,',', 1306 & f15.8,',',f15.8,')') 1307 endif 1308 endif 1309c ----- store in par_arr ------- END 1310 par(cbuf)=par(cbuf)+ 1311 & val_oo*ppm+val_ov*ppm 1312 ind=ind+4 1313 enddo ! end-loop-ix 1314 enddo ! end-loop-iy 1315 enddo ! end-loop-iatom 1316 else ! Nonrel calc 1317 if (.not. ga_create(mt_dbl,nbf,nbf, 1318 & 'gd2p1: g_t1',0,0,g_t1)) 1319 $ call errquit('get_par: g_t1',0,GA_ERR) 1320 call ga_zero(g_t1) 1321 if (.not. ga_create(mt_dbl,nbf,nbf, 1322 & 'gd2p1: g_t2',0,0,g_t2)) 1323 $ call errquit('get_par: g_t2',0,GA_ERR) 1324 call ga_zero(g_t2) 1325 if (.not. ga_create(mt_dbl,nbf,nbf, 1326 & 'gd2p1: g_t3',0,0,g_t3)) 1327 $ call errquit('get_par: g_t3',0,GA_ERR) 1328 call ga_zero(g_t3) 1329c -------- get R_nu ------------- START 1330 dims(1) =nbf 1331 chunk(1)=nbf 1332 do i=1,3 1333 if (.not. nga_create(mt_dbl,1,dims,"Array R",chunk,g_R(i))) 1334 $ call errquit('get_par: g_R', 0,GA_ERR) 1335 enddo 1336 1337 call get_chi_centers_ga(g_R,basis,nbf,geom,natoms_tot) 1338 1339! if (ga_nodeid().eq.0) 1340! & write(*,*) '----------g_R--------- START' 1341! call ga_print(g_R(1)) 1342! call ga_print(g_R(2)) 1343! call ga_print(g_R(3)) 1344! if (ga_nodeid().eq.0) 1345! & write(*,*) '----------g_R--------- END' 1346c -------- get R_nu ------------- END 1347 if (ga_nodeid().eq.0) 1348 & write(*,*) 'Calc. par tensor-> nonrel' 1349c ----- Calculate g_h01 --------- START 1350 alo(1) = nbf 1351 alo(2) = -1 1352 alo(3) = -1 1353 ahi(1) = nbf 1354 ahi(2) = nbf 1355 ahi(3) = 3*natoms_slc 1356 if (.not.nga_create(MT_DBL,3,ahi,'H01 matrix',alo,g_h01)) call 1357 & errquit('get_par: nga_create failed g_h01',0,GA_ERR) 1358 call ga_zero(g_h01) 1359 call int_giao_1ega(basis,basis,g_h01,'h01',pos(1), 1360 & natoms_slc,oskel) 1361 1362 if (debug_get_par.eq.1) then 1363 if (ga_nodeid().eq.0) 1364 & write(*,*) '------ g_h01 --------- START' 1365 call ga_print(g_h01) 1366 if (ga_nodeid().eq.0) 1367 & write(*,*) '------ g_h01 --------- END' 1368 endif 1369c ----- Calculate g_h01 --------- END 1370 alo(1) = 1 1371 ahi(1) = nbf 1372 alo(2) = 1 1373 ahi(2) = nbf 1374 blo(1) = 1 1375 bhi(1) = nbf 1376 blo(2) = 1 1377 bhi(2) = nbf 1378 blo(3) = 0 1379 bhi(3) = 0 1380 dlo(1) = 1 1381 dhi(1) = nbf 1382 dlo(2) = 1 1383 dhi(2) = nbf 1384 ind=0 1385 1386c if (ga_nodeid().eq.0) 1387c & write(*,*) '-------g_dens(1)----START' 1388c call ga_print(g_dens(1)) 1389c if (ga_nodeid().eq.0) 1390c & write(*,*) '-------g_dens(1)----END' 1391c if (ga_nodeid().eq.0) 1392c & write(*,*) '-------g_dens(2)----START' 1393c call ga_print(g_dens(2)) 1394c if (ga_nodeid().eq.0) 1395c & write(*,*) '-------g_dens(2)----END' 1396 1397 do iatom = 1, natoms_slc 1398 do iy = 1, 3 1399 blo(3) = blo(3) + 1 1400 bhi(3) = bhi(3) + 1 1401 do ix = 1, 3 1402 alo(3) = ix 1403 ahi(3) = ix 1404 val_oo = nga_ddot_patch(g_d1_oo ,'n',alo,ahi, 1405 & g_h01,'n',blo,bhi) 1406 val_ov = nga_ddot_patch(g_d1_ov ,'n',alo,ahi, 1407 & g_h01,'n',blo,bhi) 1408 cbuf=9*(iatom-1)+3*(ix-1)+iy ! transpose 1409c ----- Calc. Fukui term ------- START 1410 a=ind_tmn(1,ix) 1411 b=ind_tmn(2,ix) 1412 call nga_copy_patch('n',g_h01,blo,bhi, 1413 & g_t1 ,dlo,dhi) 1414 call ga_scale_cols(g_t1,g_R(b)) ! R_{nu,b} g_N 1415 call ga_scale_rows(g_t1,g_R(a)) ! R_{mu,a} [R_{nu,b} g_N] -> g_t1 1416 call nga_copy_patch('n',g_h01,blo,bhi, 1417 & g_t2 ,dlo,dhi) 1418 call ga_scale_cols(g_t2,g_R(a)) ! R_{nu,a} g_N 1419 call ga_scale_rows(g_t2,g_R(b)) ! R_{mu,b} [R_{nu,a} g_N] -> g_t2 1420 call ga_add(1.0d0,g_t1,-1.0d0,g_t2,g_t3) ! g_t3=-4c^2 <chi_mu|i/(2c)(R_mu x R_nu)_k h_t^{01}|chi_nu> 1421 1422c if (ga_nodeid().eq.0) then 1423c write(*,2) iatom,ix,iy 1424c 2 format('-------g_t3(',i4,',',i4,',',i4,')----START') 1425c endif 1426c call ga_print(g_t3) 1427c if (ga_nodeid().eq.0) then 1428c write(*,3) iatom,ix,iy 1429c 3 format('-------g_t3(',i4,',',i4,',',i4,')----END') 1430c endif 1431 1432 val=0.0d0 1433 do ispin=1,npol 1434 val=val+nga_ddot_patch(g_dens(ispin),'n',dlo,dhi, 1435 & g_t3,'n',dlo,dhi) 1436 enddo ! end-loop-ispin 1437 Fukui(ix,iy,iatom)=val*ppm 1438 1439c WARNING: If I attempt to print lines below the values for dia 1440c are affected somehow. This is a weird problem 1441c Otherwise, everything is fine. 1442c if (ga_nodeid().eq.0) then 1443c write(*,1) ix,iy,iatom,Fukui(ix,iy,iatom) 1444c 1 format('Fukui(',i4,',',i4,',',i4,')=',f15.8) 1445c endif 1446 1447c ----- Calc. Fukui term ------- END 1448c ----- store in par_arr ------- START 1449 par_arr(ind+1)=Fukui(ix,iy,iatom) 1450 par_arr(ind+2)=val_oo*ppm 1451 par_arr(ind+3)=val_ov*ppm 1452 par_arr(ind+4)=par_arr(ind+1)+ 1453 & val_oo*ppm+ 1454 & val_ov*ppm 1455 1456 if (debug_get_par.eq.1) then 1457 if (ga_nodeid().eq.0) then 1458 write(*,12) iatom,iy,ix, 1459 & alo(1),alo(2),alo(3), 1460 & ahi(1),ahi(2),ahi(3), 1461 & blo(1),blo(2),blo(3), 1462 & bhi(1),bhi(2),bhi(3), 1463 & par_arr(ind+1), 1464 & par_arr(ind+2),par_arr(ind+3), 1465 & par_arr(ind+4) 1466 12 format('(iatom,iy,ix)=(',i3,',',i3,',',i3,') ', 1467 & 'alo=(',i3,',',i3,',',i3,') ', 1468 & 'ahi=(',i3,',',i3,',',i3,') ', 1469 & 'blo=(',i3,',',i3,',',i3,') ', 1470 & 'bhi=(',i3,',',i3,',',i3,') ', 1471 & 'para=(',f15.8,',',f15.8,',', 1472 & f15.8,',',f15.8,')') 1473 endif 1474 endif 1475 1476c ----- store in par_arr ------- END 1477 par(cbuf)=par_arr(ind+1)+ 1478 & val_oo*ppm+val_ov*ppm 1479 ind=ind+4 1480 enddo ! end-loop-ix 1481 enddo ! end-loop-iy 1482 enddo ! end-loop-iatom 1483 if (.not.ga_destroy(g_h01)) call 1484 & errquit('get_par: ga_destroy failed g_h01',0,GA_ERR) 1485 if (.not.ga_destroy(g_t1)) call 1486 & errquit('get_par: ga_destroy failed g_t1',0,GA_ERR) 1487 if (.not.ga_destroy(g_t2)) call 1488 & errquit('get_par: ga_destroy failed g_t2',0,GA_ERR) 1489 if (.not.ga_destroy(g_t3)) call 1490 & errquit('get_par: ga_destroy failed g_t3',0,GA_ERR) 1491 do i=1,3 1492 if (.not.ga_destroy(g_R(i))) call 1493 & errquit('get_par: ga_destroy failed g_R',0,GA_ERR) 1494 enddo 1495 endif 1496c -------- symmetrize total par ------------ START 1497 do iatom = 1, natoms_slc 1498 do iy = 1, 3 1499 do ix = iy+1, 3 1500 cbuf1=9*(iatom-1)+3*(ix-1)+iy ! transpose 1501 cbuf2=9*(iatom-1)+3*(iy-1)+ix 1502 val=(par(cbuf1)+par(cbuf2))/2.0d0 1503 par(cbuf1)=val 1504 par(cbuf2)=val 1505 enddo ! end-loop-ix 1506 enddo ! end-loop-iy 1507 enddo ! end-loop-iatom 1508c -------- symmetrize total par ------------ END 1509 return 1510 end 1511c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1512c ============== get_par_JK() ===========================START 1513c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1514 subroutine get_par_JK( 1515 & par, ! OUT : paramagnetic tensor 1516 & ga_para1, ! IN : paramagnetic tensor (GA) 1517 & ga_h01_num, ! IN : H01 matrix 1518 & par_arr, ! OUT : par tensor gauge,OO,OV 1519 & g_d1_oo, ! IN : Perturbed density matrix (OO part) 1520 & g_d1_ov, ! IN : Perturbed density matrix (OV part) 1521 & g_d1_ov_Coul,! IN : Perturbed density matrix (OV part-Coul only) 1522 & g_d1_ov_Exch,! IN : Perturbed density matrix (OV part-Exch only) 1523 & g_d1_ov_noJK,! IN : Perturbed density matrix (OV part-all -(Coul,Exch)) 1524 & g_d1_ov_1e, ! IN : Perturbed density matrix (OV part-1e contrib) 1525 & g_d1_ov_eSji,! IN : Perturbed density matrix (OV part-eSji contrib) 1526 & g_dens, ! IN : e-density 1527 & npol, ! IN : nr. polarizations 1528 & Fukui, ! OUT : Fukui term 1529 & basis, ! IN : basis handle 1530 & nbf, ! IN : nr. basis functions 1531 & pos, ! IN : Nuclear positions (x,y,z) 1532 & natoms_slc, ! IN : nr. atoms selected 1533 & natoms_tot, ! IN : nr. atoms total 1534 & geom, 1535 & rtdb, 1536 & oskel) 1537c Purpose : Assemble NMR Chemical Shielding: paramagnetic tensor 1538c Author : Fredy Aquino 1539c Date : 03-03-11, 11-23-12 (adding Fukui term) 1540c Note.- Adaptation of hnd_giaox.F to hold 1541c - Nonrel(HF or DFT) in unrestricted/restricted calculations. 1542c - Relativistic (zora) in unrestricted/restricted calculations. 1543c 1544 implicit none 1545c 1546#include "errquit.fh" 1547#include "global.fh" 1548#include "mafdecls.fh" 1549#include "msgids.fh" 1550#include "rtdb.fh" 1551#include "apiP.fh" 1552#include "prop.fh" 1553#include "bgj.fh" 1554#include "bas.fh" 1555#include "stdio.fh" 1556#include "zora.fh" 1557c 1558 integer rtdb,geom 1559 integer basis 1560 integer alo(3), ahi(3), blo(3), bhi(3), 1561 & dlo(3),dhi(3) 1562 integer ld(2),cbuf,cbuf1,cbuf2 1563 integer nbf,natoms_slc,natoms_tot 1564 integer g_dens(3) ! IN : electronic density 1565 double precision Fukui(3,3,natoms_slc) ! used in analyt shield 1566 integer npol ! OUT: Used only in analytical shieldings 1567 integer i,dims(3),chunk(3), 1568 & g_R(3),ispin,a,b 1569 double precision pos(3*natoms_slc),val 1570 integer iatom,ix,iy,ixy,ind 1571 double precision val_oo,val_ov 1572 double precision par_arr(*) 1573 double precision val_ov_Coul,val_ov_Exch,val_ov_noJK, 1574 & val_ov_1e,val_ov_eSji 1575 integer g_h01, ! local 1576 & ga_h01_num 1577 integer g_d1_oo,g_d1_ov ! input: Perturbed density matrics OO and OV 1578 ! OO, due to occupied-occupied MOs 1579 ! OV, due to occupied-viritual MOs 1580 integer g_d1_ov_Coul,g_d1_ov_Exch,g_d1_ov_noJK, 1581 & g_d1_ov_1e,g_d1_ov_eSji 1582 integer ga_para1,g_t1,g_t2,g_t3 1583 integer debug_get_par 1584 double precision par(*) 1585 logical oskel 1586 integer ind_tmn(2,3) 1587 data ind_tmn / 2, 3, ! tmn=123 1588 & 3, 1, ! tmn=231 1589 & 1, 2 / ! tmn=312 1590 external get_chi_centers_ga 1591 double precision ppm 1592 data ppm /26.62566914d+00/ 1593 1594 debug_get_par=0 1595 1596 do ixy = 1, 9*natoms_slc 1597 par(ixy) = 0.0d0 ! initialize 1598 enddo 1599 if (do_zora) then 1600 if (ga_nodeid().eq.0) 1601 & write(*,*) 'Calc. par tensor-> zora' 1602c +++++ Initialize dbl_mb(k_para) with ga_para1 ++++ START 1603c ---- STORE: g_para1 --> dbl_mb(k_para) 1604 alo(1)=1 1605 ahi(1)=3 1606 alo(2)=1 1607 ahi(2)=3 1608 alo(3)=1 1609 ahi(3)=natoms_slc 1610 ld(1)=3 1611 ld(2)=3 1612 call nga_get(ga_para1,alo,ahi,par(1),ld) 1613c +++++ Initialize dbl_mb(k_para) with g_para1 ++++ END 1614 alo(1) = 1 1615 ahi(1) = nbf 1616 alo(2) = 1 1617 ahi(2) = nbf 1618 blo(1) = 1 1619 bhi(1) = nbf 1620 blo(2) = 1 1621 bhi(2) = nbf 1622 ixy = 0 1623 blo(3) = 0 1624 bhi(3) = 0 1625 ind=0 1626 do iatom = 1, natoms_slc 1627 do iy = 1, 3 1628 blo(3) = blo(3) + 1 1629 bhi(3) = bhi(3) + 1 1630 do ix = 1, 3 1631 alo(3) = ix 1632 ahi(3) = ix 1633 ixy = ixy + 1 1634c val = nga_ddot_patch(g_d1 ,'n',alo,ahi, 1635c & ga_h01_num,'n',blo,bhi) 1636 val_oo = nga_ddot_patch(g_d1_oo ,'n',alo,ahi, 1637 & ga_h01_num,'n',blo,bhi) 1638 val_ov = nga_ddot_patch(g_d1_ov ,'n',alo,ahi, 1639 & ga_h01_num,'n',blo,bhi) 1640 val_ov_Coul = nga_ddot_patch(g_d1_ov_Coul,'n',alo,ahi, 1641 & ga_h01_num,'n',blo,bhi) 1642 val_ov_Exch = nga_ddot_patch(g_d1_ov_Exch,'n',alo,ahi, 1643 & ga_h01_num,'n',blo,bhi) 1644 val_ov_noJK = nga_ddot_patch(g_d1_ov_noJK,'n',alo,ahi, 1645 & ga_h01_num,'n',blo,bhi) 1646 val_ov_1e = nga_ddot_patch(g_d1_ov_1e,'n',alo,ahi, 1647 & ga_h01_num,'n',blo,bhi) 1648 val_ov_eSji = nga_ddot_patch(g_d1_ov_eSji,'n',alo,ahi, 1649 & ga_h01_num,'n',blo,bhi) 1650 1651 cbuf=9*(iatom-1)+3*(ix-1)+iy !-1 transpose 1652c ----- store in par_arr ------- START 1653 par_arr(ind+1)=par(cbuf) 1654 par_arr(ind+2)=val_oo*ppm 1655 par_arr(ind+3)=val_ov*ppm 1656 par_arr(ind+4)=par(cbuf)+ 1657 & val_oo*ppm+ 1658 & val_ov*ppm 1659 par_arr(ind+5)=val_ov_Coul*ppm 1660 par_arr(ind+6)=val_ov_Exch*ppm 1661 par_arr(ind+7)=val_ov_noJK*ppm 1662 par_arr(ind+8)=val_ov_1e*ppm 1663 par_arr(ind+9)=val_ov_eSji*ppm 1664c ----- store in par_arr ------- END 1665 par(cbuf)=par(cbuf)+ 1666 & val_oo*ppm+val_ov*ppm 1667 ind=ind+9 1668 enddo ! end-loop-ix 1669 enddo ! end-loop-iy 1670 enddo ! end-loop-iatom 1671 else ! Nonrel calc 1672 if (.not. ga_create(mt_dbl,nbf,nbf, 1673 & 'gd2p1: g_t1',0,0,g_t1)) 1674 $ call errquit('get_par_JK: g_t1',0,GA_ERR) 1675 call ga_zero(g_t1) 1676 if (.not. ga_create(mt_dbl,nbf,nbf, 1677 & 'gd2p1: g_t2',0,0,g_t2)) 1678 $ call errquit('get_par_JK: g_t2',0,GA_ERR) 1679 call ga_zero(g_t2) 1680 if (.not. ga_create(mt_dbl,nbf,nbf, 1681 & 'gd2p1: g_t3',0,0,g_t3)) 1682 $ call errquit('get_par_JK: g_t3',0,GA_ERR) 1683 call ga_zero(g_t3) 1684c -------- get R_nu ------------- START 1685 dims(1) =nbf 1686 chunk(1)=nbf 1687 do i=1,3 1688 if (.not. nga_create(mt_dbl,1,dims,"Array R",chunk,g_R(i))) 1689 $ call errquit('get_par_JK: g_R', 0,GA_ERR) 1690 enddo 1691 call get_chi_centers_ga(g_R,basis,nbf,geom,natoms_tot) 1692c -------- get R_nu ------------- END 1693 if (ga_nodeid().eq.0) 1694 & write(*,*) 'Calc. par tensor-> nonrel' 1695c ----- Calculate g_h01 --------- START 1696 alo(1) = nbf 1697 alo(2) = -1 1698 alo(3) = -1 1699 ahi(1) = nbf 1700 ahi(2) = nbf 1701 ahi(3) = 3*natoms_slc 1702 if (.not.nga_create(MT_DBL,3,ahi,'H01 matrix',alo,g_h01)) call 1703 & errquit('get_par_JK: nga_create failed g_h01',0,GA_ERR) 1704 call ga_zero(g_h01) 1705 call int_giao_1ega(basis,basis,g_h01,'h01',pos(1), 1706 & natoms_slc,oskel) 1707 if (debug_get_par.eq.1) then 1708 if (ga_nodeid().eq.0) 1709 & write(*,*) '------ g_h01 --------- START' 1710 call ga_print(g_h01) 1711 if (ga_nodeid().eq.0) 1712 & write(*,*) '------ g_h01 --------- END' 1713 endif 1714c ----- Calculate g_h01 --------- END 1715 alo(1) = 1 1716 ahi(1) = nbf 1717 alo(2) = 1 1718 ahi(2) = nbf 1719 blo(1) = 1 1720 bhi(1) = nbf 1721 blo(2) = 1 1722 bhi(2) = nbf 1723 blo(3) = 0 1724 bhi(3) = 0 1725 ind=0 1726 do iatom = 1, natoms_slc 1727 do iy = 1, 3 1728 blo(3) = blo(3) + 1 1729 bhi(3) = bhi(3) + 1 1730 do ix = 1, 3 1731 alo(3) = ix 1732 ahi(3) = ix 1733 val_oo = nga_ddot_patch(g_d1_oo ,'n',alo,ahi, 1734 & g_h01,'n',blo,bhi) 1735 val_ov = nga_ddot_patch(g_d1_ov ,'n',alo,ahi, 1736 & g_h01,'n',blo,bhi) 1737 cbuf=9*(iatom-1)+3*(ix-1)+iy ! transpose 1738c ----- Calc. Fukui term ------- START 1739 a=ind_tmn(1,ix) 1740 b=ind_tmn(2,ix) 1741 call nga_copy_patch('n',g_h01,blo,bhi, 1742 & g_t1 ,dlo,dhi) 1743 call ga_scale_cols(g_t1,g_R(b)) ! R_{nu,b} g_N 1744 call ga_scale_rows(g_t1,g_R(a)) ! R_{mu,a} [R_{nu,b} g_N] -> g_t1 1745 call nga_copy_patch('n',g_h01,blo,bhi, 1746 & g_t2 ,dlo,dhi) 1747 call ga_scale_cols(g_t2,g_R(a)) ! R_{nu,b} g_N 1748 call ga_scale_rows(g_t2,g_R(b)) ! R_{mu,a} [R_{nu,b} g_N] -> g_t2 1749 call ga_add(1.0d0,g_t1,-1.0d0,g_t2,g_t3) ! g_t3=-4c^2 <chi_mu|i/(2c)(R_mu x R_nu)_k h_t^{01}|chi_nu> 1750 val=0.0d0 1751 do ispin=1,npol 1752 val=val+nga_ddot_patch(g_dens(ispin),'n',dlo,dhi, 1753 & g_t3,'n',dlo,dhi) 1754 enddo ! end-loop-ispin 1755 1756 Fukui(ix,iy,iatom)=val*ppm 1757 1758c write(*,1) ix,iy,iatom,Fukui(ix,iy,iatom) 1759c 1 format('Fukui(',i4,',',i4,',',i4,')=',f15.8) 1760 1761c ----- Calc. Fukui term ------- END 1762c ----- store in par_arr ------- START 1763 par_arr(ind+1)=Fukui(ix,iy,iatom) 1764 par_arr(ind+2)=val_oo*ppm 1765 par_arr(ind+3)=val_ov*ppm 1766 par_arr(ind+4)=par_arr(ind+1)+ 1767 & val_oo*ppm+ 1768 & val_ov*ppm 1769 if (debug_get_par.eq.1) then 1770 if (ga_nodeid().eq.0) then 1771 write(*,12) iatom,iy,ix, 1772 & alo(1),alo(2),alo(3), 1773 & ahi(1),ahi(2),ahi(3), 1774 & blo(1),blo(2),blo(3), 1775 & bhi(1),bhi(2),bhi(3), 1776 & par_arr(ind+1), 1777 & par_arr(ind+2),par_arr(ind+3), 1778 & par_arr(ind+4) 1779 12 format('(iatom,iy,ix)=(',i3,',',i3,',',i3,') ', 1780 & 'alo=(',i3,',',i3,',',i3,') ', 1781 & 'ahi=(',i3,',',i3,',',i3,') ', 1782 & 'blo=(',i3,',',i3,',',i3,') ', 1783 & 'bhi=(',i3,',',i3,',',i3,') ', 1784 & 'para=(',f15.8,',',f15.8,',', 1785 & f15.8,',',f15.8,')') 1786 endif 1787 endif 1788c ----- store in par_arr ------- END 1789 par(cbuf)=par_arr(ind+1)+ 1790 & val_oo*ppm+val_ov*ppm 1791 ind=ind+4 1792 enddo ! end-loop-ix 1793 enddo ! end-loop-iy 1794 enddo ! end-loop-iatom 1795 if (.not.ga_destroy(g_h01)) call 1796 & errquit('get_par_JK: ga_destroy failed g_h01',0,GA_ERR) 1797 if (.not.ga_destroy(g_t1)) call 1798 & errquit('get_par_JK: ga_destroy failed g_t1',0,GA_ERR) 1799 if (.not.ga_destroy(g_t2)) call 1800 & errquit('get_par_JK: ga_destroy failed g_t2',0,GA_ERR) 1801 if (.not.ga_destroy(g_t3)) call 1802 & errquit('get_par_JK: ga_destroy failed g_t3',0,GA_ERR) 1803 do i=1,3 1804 if (.not.ga_destroy(g_R(i))) call 1805 & errquit('get_par_JK: ga_destroy failed g_R',0,GA_ERR) 1806 enddo 1807 endif 1808c -------- symmetrize total par ------------ START 1809 do iatom = 1, natoms_slc 1810 do iy = 1, 3 1811 do ix = iy+1, 3 1812 cbuf1=9*(iatom-1)+3*(ix-1)+iy ! transpose 1813 cbuf2=9*(iatom-1)+3*(iy-1)+ix 1814 val=(par(cbuf1)+par(cbuf2))/2.0d0 1815 par(cbuf1)=val 1816 par(cbuf2)=val 1817 enddo ! end-loop-ix 1818 enddo ! end-loop-iy 1819 enddo ! end-loop-iatom 1820c -------- symmetrize total par ------------ END 1821 return 1822 end 1823c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1824c ============== get_par_JK() =============================END 1825c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1826 subroutine get_dia(dia, ! OUT : paramagnetic tensor 1827 & ga_dia, ! IN : dia tensor 1828 & basis, ! IN : basis handle 1829 & g_dens, ! IN : e-density 1830 & nbf, ! IN : nr. basis functions 1831 & npol, ! IN : nr. of polarizations 1832 & Fukui, ! IN : Fukui term 1833 & pos, ! IN : Nuclear positions (x,y,z) 1834 & natoms, ! IN : nr. atoms 1835 & rtdb, 1836 & oskel) 1837 1838c Purpose : Assemble NMR Chemical Shielding: diamagnetic tensor 1839c Author : Fredy Aquino 1840c Date : 03-03-11, 11-23-12 (adding Fukui term) 1841c Note.- Adaptation of hnd_giaox.F to hold 1842c - Nonrel(HF or DFT) in unrestricted/restricted calculations. 1843c - Relativistic (zora) in unrestricted/restricted calculations. 1844c 1845 implicit none 1846c 1847#include "errquit.fh" 1848#include "global.fh" 1849#include "mafdecls.fh" 1850#include "msgids.fh" 1851#include "rtdb.fh" 1852#include "bas.fh" 1853#include "apiP.fh" 1854#include "prop.fh" 1855#include "bgj.fh" 1856#include "stdio.fh" 1857#include "zora.fh" 1858c 1859c Global arrays variables defined in zora.fh 1860c --> ga_dia,ga_para1,ga_Fji,ga_h01_num 1861 integer rtdb 1862 integer basis ! [input] basis handle 1863 integer g_dens(3) ! IN: electronic density 1864 integer ga_dia 1865 integer alo(3), ahi(3), blo(3), bhi(3) 1866 integer ld(2),cbuf, 1867 & ix1,iy1 1868 integer nbf,natoms,npol,ispin 1869 double precision ac,pos(3*natoms) 1870 double precision Fukui(3,3,natoms) ! used in analyt shield 1871 integer iatom,ix,iy,ixy,ind1,ind2 1872 integer g_h11 ! local ga 1873 double precision dia(*) 1874 logical oskel 1875 double precision ppm,val 1876 data ppm /26.62566914d+00/ 1877 do ixy = 1, 9*natoms 1878 dia(ixy) = 0.0d0 ! initialize 1879 enddo 1880 if (do_zora) then ! get-dia---- START 1881c ---- STORE: ga_dia --> dbl_mb(k_dia) 1882 alo(1)=1 1883 ahi(1)=3 1884 alo(2)=1 1885 ahi(2)=3 1886 alo(3)=1 1887 ahi(3)=natoms 1888 ld(1)=3 1889 ld(2)=3 1890 call nga_get(ga_dia,alo,ahi,dia(1),ld) 1891 else 1892 alo(1) = nbf 1893 alo(2) = -1 1894 alo(3) = -1 1895 ahi(1) = nbf 1896 ahi(2) = nbf 1897 ahi(3) = 9*natoms 1898 if (.not.nga_create(MT_DBL,3,ahi,'H11 matrix',alo,g_h11)) call 1899 & errquit('get_dia: nga_create failed g_h11 all',0,GA_ERR) 1900 call ga_zero(g_h11) 1901 call int_giao_1ega(basis,basis,g_h11,'h11 para', 1902 & pos(1),natoms,oskel) 1903 blo(1) = 1 1904 bhi(1) = nbf 1905 blo(2) = 1 1906 bhi(2) = nbf 1907 blo(3) = 0 1908 bhi(3) = 0 1909 alo(1) = 1 1910 ahi(1) = nbf 1911 alo(2) = 1 1912 ahi(2) = nbf 1913 ixy = 0 1914 do iatom = 1, natoms 1915 ix1=1 1916 iy1=1 1917 do ix = 1, 9 1918 if (ix1.gt.3) then 1919 ix1=1 1920 iy1=iy1+1 1921 endif 1922 ixy = ixy + 1 1923 alo(3) = ixy 1924 ahi(3) = ixy 1925 do ispin=1,npol 1926 val=nga_ddot_patch(g_dens(ispin),'n',blo,bhi, 1927 & g_h11,'n',alo,ahi) 1928 dia(ixy)=dia(ixy) + val * ppm 1929 enddo ! end-loop-ispin 1930 dia(ixy)=dia(ixy)-Fukui(ix1,iy1,iatom) 1931c Note.- Printing lines below does not affect output 1932c It is ok to print it out but I will comment to 1933c reduce printouts. 1934 1935c if (ga_nodeid().eq.0) then 1936c write(*,1) ix1,iy1,iatom,Fukui(ix1,iy1,iatom) 1937c 1 format('Fukui(ix1,iy1,iatom)(',i4,',',i4,',',i4,')=', 1938c & f15.8) 1939c endif 1940 1941 ix1=ix1+1 1942 enddo ! end-loop-ix 1943 enddo ! end-loop-atoms 1944c 1945c s(dia)xy = s(dia)xy + Sum(n,l) D0(n,l) * H11(dia)xy(n,l) 1946 1947 call ga_zero(g_h11) 1948 call int_giao_1ega(basis,basis,g_h11,'h11 dia', 1949 & pos(1),natoms,oskel) 1950 alo(1) = 1 1951 ahi(1) = nbf 1952 alo(2) = 1 1953 ahi(2) = nbf 1954 ixy=0 1955 do iatom = 1, natoms 1956 ix1=1 1957 iy1=1 1958 do ix = 1, 9 1959 if (ix1.gt.3) then 1960 ix1=1 1961 iy1=iy1+1 1962 endif 1963 ixy = ixy + 1 1964 alo(3) = ixy 1965 ahi(3) = ixy 1966 ac=0.0d0 1967 do ispin=1,npol 1968 val=nga_ddot_patch(g_dens(ispin),'n',blo,bhi, 1969 & g_h11,'n',alo,ahi) 1970 ac=ac + val * ppm 1971 enddo ! end-loop-ispin 1972 dia(ixy)=dia(ixy)+ac 1973c 1974c WARNING: If I try to print ac, uncommenting lines below 1975c the dia(ixy) values are affected somehow. (FA-11-24-12) 1976c This is a weird problem, otherwise everything comes 1977c fine. 1978c if (ga_nodeid().eq.0) then 1979c write(*,2) ix1,iy1,iatom,ac 1980c 2 format('h11-dia(',i4,',',i4,',',i4,')=', 1981c & f15.8) 1982c endif 1983 1984 ix1=ix1+1 1985 enddo 1986 enddo 1987 if (.not.ga_destroy(g_h11)) call 1988 & errquit('get_dia: ga_destroy failed g_h11',0,GA_ERR) 1989 endif ! get-dia---- END 1990c ---------- symmetrize dia ----------- START 1991 do iatom=1,natoms 1992 do ix=1,3 1993 do iy=ix+1,3 1994 ind1=9*(iatom-1)+3*(iy-1)+ix 1995 ind2=9*(iatom-1)+3*(ix-1)+iy 1996 val=(dia(ind1)+dia(ind2))/2.0d0 1997 dia(ind1)=val 1998 dia(ind2)=val 1999 enddo ! end-loop-iy 2000 enddo ! end-loop ix 2001 enddo ! end-loop-iat 2002c ---------- symmetrize dia ----------- END 2003 return 2004 end 2005c -------------- for shieldings NMLO analysis ----------------- START 2006 subroutine create_munu4nbo_shield( 2007 & rtdb, ! in: rtdb handle 2008 & g_tvec, ! in: eigenvectors or T diagonalizing matrix 2009 & nat, ! in: nr. atoms 2010 & atmnr, ! in: selected atoms 2011 & basis, ! in: basis handle 2012 & npol, ! in: nr. polarizaitons 2013 & nocc, ! in: nr. occ nocc(i) i=1,npol 2014 & nvirt, ! in: nr. virt nvirt(i) i=1,npol 2015 & nmo) ! in: nr. MO 2016c 2017 implicit none 2018c 2019#include "nwc_const.fh" 2020#include "errquit.fh" 2021#include "global.fh" 2022#include "bas.fh" 2023#include "mafdecls.fh" 2024#include "geom.fh" 2025#include "stdio.fh" 2026#include "rtdb.fh" 2027#include "cosmo.fh" 2028#include "msgids.fh" 2029#include "zora.fh" 2030c 2031c FA: Revised on 06-22-11 2032c ------Main outputs -------- START 2033 integer g_munu_rot, ! hyp-FCSD dia part 2034 & g_munu_rot1, 2035 & g_munu_rot2, 2036 & g_acc2 ! hyp-PSOSO para part 2037c ------Main outputs -------- END 2038 integer rtdb,basis 2039 integer g_dens,g_tvec,atmnr(nat) 2040 integer g_munuCSdia,g_munuCSpar1,g_munuCSHpar, 2041 & g_tnp,g_acc,vectors(2), 2042 & g_tnp1,g_acc1,g_acc3 2043 integer vectors_scl(2),ispin,iat1 2044 double precision ac_val,val1,sign 2045 integer i,j,k,m,n,ndir,ndir1 2046 integer jlo,jhi,s,nbf,nmo,nsize,nsize1,nsize2 2047 integer npol,npol_munu,ntot 2048 integer ind,nlst,count,nocc(npol),nvirt(npol) 2049 integer Natoms_munu,atmnr_munu(nat) 2050 integer Ndir_munu 2051c Ndir_munu, Nr. of directions stored 2052c =3 xx yy zz 2053 double precision coeff,fact,para_rot(9) 2054 double precision tmn(2),chcdata(3) 2055 integer jlo1,jhi1,jlo2,jhi2 2056 integer g_dens1,g_c1, 2057 & g_p10,g_munuCSHpar2d,g_munuCSpar12d 2058 integer iind(2),jind(2),icalczora,type_NMR,iat,nat 2059 integer alo(3),ahi(3),elo(3),ehi(3),flo(3),fhi(3) 2060 logical dft_zoraShield_NLMOAnalysis_read ! for read-nlmo-mat 2061 character*255 zorafilename ! for read-nlmo-mat 2062 integer arr_ind(6,2) 2063 data ((arr_ind(j,i),i=1,2),j=1,6) 2064 & /1,1,2,2,3,3,1,2,1,3,2,3/ 2065 external dft_zoraShield_NLMOAnalysis_read,get_P10_rot, 2066 & fill_munuPSOSO_1,get_par_gshift_rot, 2067 & get2dmat,wshldfile 2068c --> To store ONLY munu principal components xx,yy,zz 2069c g_munuCSdia is created in dft_zora_EPR.F 2070c size(g_munuCSdia)=nlst*ndir (linear vector) 2071c g_dens, spin density matrix 2072c nbf x nbf elements (bidimensional matrix) 2073c Legend: 2074c ndir=6 = xx, yy, zz, xy, xz, yz 2075c nbf, Nr of basis functions 2076c nlst=nbf*(nbf+1)/2 2077 2078 if (.not. bas_numbf(basis,nbf)) call errquit 2079 & ('munu: bas_numbf failed',555, BASIS_ERR) 2080 Natoms_munu=nat 2081 do i=1,Natoms_munu 2082 atmnr_munu(i)=atmnr(i) 2083 enddo 2084 npol_munu=npol 2085 Ndir_munu=3 2086 nlst=nbf*(nbf+1)/2 ! size of xx,yy,zz,xy,xz,yz chunk 2087c ++++++ Read NLMO matrices +++++++++ START 2088 ndir =6 2089 ndir1=3 2090 call util_file_name(lbl_nlmoshield,.false.,.false., 2091 & zorafilename) 2092 icalczora = 0 ! initialize the flag 2093 if (.not.dft_zoraShield_NLMOAnalysis_read( 2094 & zorafilename, ! in : filename 2095 & nbf, ! in : nr basis functions 2096 & ndir, ! in : nr of directions: 6 = xx yy zz xy xz yz 2097 & ndir1, ! in : nr of directions: 3 = x y z 2098 & nat, ! in : list of selected atoms 2099 & nocc, ! in : nocc(i) i=1,2 nr. occupations in alpha and beta 2100 & npol, ! in: nr polarizations 2101 & g_munuCSdia, ! out: munu matrix of dia 2102 & g_munuCSpar1, ! out: munu matrix of par 1st term 2103 & g_munuCSHpar, ! out: munu matrix of H10 2104 & vectors, ! out: MOs 2105 & g_c1, ! out: perturbed MO 2106 & g_dens)) icalczora=1 2107 2108 goto 10 2109 if (ga_nodeid().eq.0) 2110 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuCSdia -- START' 2111 call ga_print(g_munuCSdia) 2112 if (ga_nodeid().eq.0) 2113 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuCSdia -- END' 2114 if (ga_nodeid().eq.0) 2115 & write(*,*) 2116 & 'AFT-reading-HYP-NLMO matrices: g_munuCSHpar -- START' 2117 call ga_print(g_munuCSHpar) 2118 if (ga_nodeid().eq.0) 2119 & write(*,*) 2120 & 'AFT-reading-HYP-NLMO matrices: g_munuCSHpar -- END' 2121 if (ga_nodeid().eq.0) 2122 & write(*,*) 2123 & 'AFT-reading-HYP-NLMO matrices: g_munuCSpar1 -- START' 2124 call ga_print(g_munuCSpar1) 2125 if (ga_nodeid().eq.0) 2126 & write(*,*) 2127 & 'AFT-reading-HYP-NLMO matrices: g_munuCSpar1 -- END' 2128 if (ga_nodeid().eq.0) 2129 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_dens ----- START' 2130 call ga_print(g_dens) 2131 if (ga_nodeid().eq.0) 2132 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_dens ----- END' 2133 if (ga_nodeid().eq.0) 2134 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_c1 ----- START' 2135 call ga_print(g_c1) 2136 if (ga_nodeid().eq.0) 2137 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_c1 ----- END' 2138 10 continue 2139c ++++++ Read NLMO matrices +++++++++ END 2140 call get_unique_elmat(g_dens,g_dens1,nlst,nbf) ! out: g_dens1 2141 ndir =6 ! Nr. of directions: xx,yy,zz,xy,xz,yz 2142 ndir1=3 ! Nr. of directions: x,y,z 2143 nsize =nbf*(nbf+1)/2 ! size of xx,yy,zz,xy,xz,yz chunk 2144 nsize1=nsize*ndir ! size of whole munu per atom 2145 nsize2=nsize*ndir1 ! size of whole munu per atom 2146 if (.not. ga_create(mt_dbl,1,nsize, 2147 & 'munu4nbo: g_tnp',0,0,g_tnp)) 2148 $ call errquit('munu4nbo: g_tnp', 0,GA_ERR) 2149 call ga_zero(g_tnp) 2150 if (.not. ga_create(mt_dbl,1,nsize, 2151 & 'munu4nbo: g_tnp1',0,0,g_tnp1)) 2152 $ call errquit('munu4nbo: g_tnp1', 0,GA_ERR) 2153 call ga_zero(g_tnp1) 2154 if (.not. ga_create(mt_dbl,1,nsize, 2155 & 'munu4nbo: g_acc',0,0,g_acc)) 2156 $ call errquit('munu4nbo: g_acc', 0,GA_ERR) 2157 call ga_zero(g_acc) 2158 2159 if (.not. ga_create(mt_dbl,1,nsize, 2160 & 'munu4nbo: g_acc3',0,0,g_acc3)) 2161 $ call errquit('munu4nbo: g_acc3', 0,GA_ERR) 2162 call ga_zero(g_acc3) 2163 if (.not. ga_create(mt_dbl,1,nsize, 2164 & 'munu4nbo: g_acc1',0,0,g_acc1)) 2165 $ call errquit('munu4nbo: g_acc1', 0,GA_ERR) 2166 call ga_zero(g_acc1) 2167 alo(1) = nbf 2168 alo(2) = -1 2169 alo(3) = -1 2170 ahi(1) = nbf 2171 ntot=nocc(1)+nocc(2) 2172 ahi(2) = ntot 2173 ahi(3) = 3*nat 2174 2175 if (.not.nga_create(MT_DBL,3,ahi,'g_acc2 matrix', 2176 & alo,g_acc2)) call 2177 & errquit('g_acc2: nga_create failed g_acc2',0,GA_ERR) 2178 call ga_zero(g_acc2) 2179 if (.not. ga_create(mt_dbl,1,nlst*3*nat, 2180 & 'munu4nbo: g_munu_rot',0,0,g_munu_rot)) 2181 $ call errquit('munu4nbo: g_munu_rot', 0,GA_ERR) 2182 if (.not. ga_create(mt_dbl,1,nlst*3*nat, 2183 & 'munu4nbo: g_munu_rot2',0,0,g_munu_rot2)) 2184 $ call errquit('munu4nbo: g_munu_rot2', 0,GA_ERR) 2185 if (.not. ga_create(mt_dbl,1,nlst*3*nat, 2186 & 'munu4nbo: g_munu_rot1',0,0,g_munu_rot1)) 2187 $ call errquit('munu4nbo: g_munu_rot1', 0,GA_ERR) 2188 alo(1) = nbf 2189 alo(2) = -1 2190 alo(3) = -1 2191 ahi(1) = nbf 2192 ahi(2) = nbf 2193 ahi(3) = 3 2194 if (.not.nga_create(mt_dbl,3,ahi,'g_munuCSHpar2d matrix', 2195 & alo,g_munuCSHpar2d)) call 2196 & errquit('munu4nbo: nga_create failed g_munuCSHpar2d', 2197 & 0,GA_ERR) 2198 call ga_zero(g_munuCSHpar2d) 2199 if (.not.nga_create(mt_dbl,3,ahi,'g_munuCSpar12d matrix', 2200 & alo,g_munuCSpar12d)) call 2201 & errquit('munu4nbo: nga_create failed g_munuCSpar12d', 2202 & 0,GA_ERR) 2203 call ga_zero(g_munuCSpar12d) 2204 if (ga_nodeid().eq.0) 2205 & write(*,*) 'CHCooooooooooooo', 2206 & ' NW-Shieldings: Summary C+HC data [ppm] ', 2207 & 'ooooooooooooooo START' 2208 do iat1=1,nat 2209 call ga_zero(g_munuCSpar12d) 2210 iat=atmnr(iat1) 2211 do n=1,3 ! xx,yy,zz 2212 m=n ! For principal components ONLY 2213c ----- Do: A'= T^t A T, calculate only [A']_pp --> (do n=1,3 m=n) 2214c a_pp'= \sum_i t_ip a_ii t_ip + 2215c 2 \sum_{j=2}^n \sum_{i=1}^{j-1} t_jp a_ji t_ip 2216c g_munu_rot = A' 2217c WARNING: g_munu_rot, contains several rotated matrices 2218c since the matrices are symmetric I store only 2219c the main diagonal + lower (upper) triangular 2220c matrix in a format that looks like: 2221c a_11 a_22 ... a_nn 2222c a_21 2223c a_31 a_32 2224c a_41 a_42 a_43 2225c ... 2226c a_n1 a_n2 ... a_{n(n-1)} 2227c There are two additional transformations on g_munu_rot 2228c before leaving this routine and entering wefgfile() 2229c 1. I make the diagonalized matrix traceless 2230c ===== Transform xx_munu to 2xx_munu-(yy_munu+zz_munu) = START 2231c or 3xx_munu-(xx_munu+yy_munu+zz_munu) 2232c 2. I need to do a reordering of elements so that it is 2233c compatible with wefgfile() 2234c call reorder_munu(g_munu_rot,nat,nlst,nbf,Ndir_munu) 2235c -------------------------------------------------------------- 2236 call ga_zero(g_acc) 2237 call ga_zero(g_acc3) 2238 do s=1,6 2239c ------- get coeff() --- START 2240 iind(1)=1 2241 iind(2)=1 2242 jind(1)=9*(iat1-1)+3*(arr_ind(s,1)-1)+m 2243 jind(2)=9*(iat1-1)+3*(arr_ind(s,2)-1)+n 2244 call ga_gather(g_tvec,tmn,iind,jind,2) 2245 fact=1.0d0 2246 if (s.gt.3) fact=2.0d0 2247 coeff=fact*tmn(1)*tmn(2) 2248c ------- get coeff() --- END 2249c Note.- g_munuFCSD will be the (hyp-diag)_uv matrix 2250 jlo=nsize1*(iat1-1)+nsize*(s-1)+1 2251 jhi=jlo+nsize-1 2252 call ga_copy_patch('n',g_munuCSdia,1,1,jlo,jhi, 2253 & g_tnp ,1,1,1 ,nsize) 2254 call ga_add(1.0d0,g_acc,coeff,g_tnp,g_acc) 2255 2256 call ga_copy_patch('n',g_munuCSpar1,1,1,jlo,jhi, 2257 & g_tnp ,1,1,1 ,nsize) 2258 call ga_add(1.0d0,g_acc3,coeff,g_tnp,g_acc3) 2259 enddo ! end-loop-s 2260 2261c Note.- g_acc = \widetilde{H}_{mu nu}^{(m,m)} 2262c it is the rotated munu matrix using: T^t H T 2263 call ga_zero(g_acc1) 2264 do s=1,3 2265c ------- get coeff() --- START 2266 iind(1)=1 2267 jind(1)=9*(iat1-1)+3*(s-1)+m 2268 call ga_gather(g_tvec,tmn,iind,jind,1) 2269c if (ga_nodeid().eq.0) then 2270c write(*,1) m,s,tmn(1) 2271c 1 format('(m,s,tvec)=(',i5,',',i5,',',f15.8,')') 2272c endif 2273 coeff=tmn(1) 2274c ------- get coeff() --- END 2275c Note.- g_munuCSHpar will be the (g-shift-para)_uv matrix 2276 jlo=nsize2*(iat1-1)+nsize*(s-1)+1 2277 jhi=jlo+nsize-1 2278 call ga_copy_patch('n',g_munuCSHpar,1,1,jlo,jhi, 2279 & g_tnp1 ,1,1,1 ,nsize) 2280 call ga_add(1.0d0,g_acc1,coeff,g_tnp1,g_acc1) 2281c ----- Calculate rotated perturbed MO: g_acc2 ----- START 2282c \sum_{s=1,3} t_{sm} C_{ri}^{(s) sigma} --> g_acc2 2283 elo(1) = 1 2284 ehi(1) = nbf 2285 elo(2) = 1 2286 ehi(2) = ntot 2287 elo(3) = s 2288 ehi(3) = s 2289 flo(1) = 1 2290 fhi(1) = nbf 2291 flo(2) = 1 2292 fhi(2) = ntot 2293 flo(3) = 3*(iat1-1)+m 2294 fhi(3) = 3*(iat1-1)+m 2295 call nga_add_patch(1.0d0,g_acc2,flo,fhi, 2296 & coeff,g_c1 ,elo,ehi, 2297 & g_acc2,flo,fhi) 2298c ----- Calculate rotated perturbed MO: g_acc2 ----- END 2299 enddo ! end-loop-s 2300c Note: g_acc1 = \widetilde{H}_{mu nu}^{(m)} m=x,y,z 2301c it is the rotated munu matrix using: T H 2302c ====== Store final munu matrices === START 2303 jlo2=nlst*Ndir_munu*(iat1-1)+ 2304 & nlst*(n-1)+1 2305 jhi2=jlo2+nlst-1 2306 call ga_copy_patch('n',g_acc ,1,1, 1,nlst, 2307 & g_munu_rot,1,1,jlo2,jhi2) 2308 call ga_copy_patch('n',g_acc3 ,1,1, 1,nlst, 2309 & g_munu_rot2,1,1,jlo2,jhi2) 2310 call ga_copy_patch('n',g_acc1 ,1,1, 1,nlst, 2311 & g_munu_rot1,1,1,jlo2,jhi2) 2312c ====== Store final munu matrices === END 2313c ++++++++++++++++++CHECK++++ DIAGONALIZATION ==== START 2314c ==== sum (g_acc .* g_dens1 + Nuclear CONTRIB) 2315c = TOTAL Shieldings diagonalized 2316 jlo1=1+nbf 2317 jhi1=nsize 2318 call ga_scale_patch(g_acc,1,1,jlo1,jhi1,2.0d0) 2319 chcdata(m)=ga_ddot(g_acc,g_dens1) 2320c ++++++++++++++++++CHECK++++ DIAGONALIZATION ==== END 2321c +++++++++++++++++++++++++++++++++++++++++++++++++++++ 2322 enddo ! end-loop-n 2323 if (ga_nodeid().eq.0) then 2324 write(*,23) iat, 2325 & chcdata(1), ! dia-x 2326 & chcdata(2),chcdata(3) ! dia-y,z 2327 23 format(' CHC dia(xx,yy,zz)(',i3,')=(', 2328 & f15.8,',',f15.8,',',f15.8,')') 2329 endif 2330c ++++++++++ CHECK diagonalization in para hyperfine ++++++++ START 2331c Note.- Variables defined in zora.fh: 2332c g_CiFull, zora scaling factors 2333c filled out with values in dft_zora_scale 2334c zora switches: do_zora,do_NonRel,not_zora_scale 2335 type_NMR=1 ! =1,2,3=shieldings,hyperfine,gshift 2336 call get_P10_rot( 2337 & g_p10, ! out: Perturbed density matrix (munu nbf x nbf x 3 square matrix) 2338 & type_NMR, ! in: =1,2,3=shieldings,hyperfine,gshift 2339 & g_acc2, ! in: rotated perturbed MO vector 2340 & vectors,g_CiFull, ! in: to build zora scaled MO vector 2341 & iat1, ! in: index for selected atom nr =1,nat 2342 & nbf,nmo,npol,nocc,nvirt, 2343 & do_zora,do_NonRel,not_zora_scale,rtdb) 2344 call fill_munuPSOSO_1( ! g_munuCSHpar --> g_munuCSHpar2d 2345 & g_munu_rot1, ! in: array with unique elements 2346 & g_munuCSHpar2d, !out: nbf x nbf x 3 munu matrix for ith atom 2347 & iat1, ! in: = 1,2,...,nat 2348 & 2, ! in: type_symm = 1 symm = 2 antisymm 2349 & nbf) 2350 call fill_munuPSOSO_1( ! g_munuCSpar1 --> g_munuCSHpar2d 2351 & g_munu_rot2, ! in: array with unique elements 2352 & g_munuCSpar12d, !out: nbf x nbf x 3 munu matrix for ith atom 2353 & iat1, ! in: = 1,2,...,nat 2354 & 1, ! in: type_symm = 1 symm = 2 antisymm 2355 & nbf) 2356c +++++++++ NOW: do ddot product to get diagonalized tensor +++ START 2357 call get_par_gshift_rot( 2358 & g_dens, ! IN : spin-density 2359 & g_munuCSpar12d, ! IN : par 1st term 2360 & g_munuCSHpar2d, ! IN : h01 matrix 2361 & g_p10, ! IN : Perturbed density matrix 2362 & basis,nbf,iat,rtdb) 2363c +++++++++ NOW: do ddot product to get diagonalized tensor +++ END 2364c ++++++++++ CHECK diagonalization in para hyperfine ++++++++ END 2365 if (.not. ga_destroy(g_p10)) call errquit( 2366 & 'create_munu4nbo_shield: ga_destroy failed g_p10',0, GA_ERR) 2367 enddo ! end-loop-iat 2368 if (ga_nodeid().eq.0) 2369 & write(*,*) 'CHCooooooooooooo', 2370 & ' NW-Shieldings: Summary C+HC data [ppm] ', 2371 & 'ooooooooooooooo END' 2372c --> Main outputs: g_acc2 ,rotated perturbed MO nbf*ntot*ndir*nat 2373c ndir=1,2,3=x,y,z 2374c nbf, nr of basis functions 2375c ntot=nocc(1)+nocc(2) 2376c nat, nr of selected atoms 2377c g_munu_rot1,rotated perturbed AO matrix 2378c storing only diag + off-diag elements 2379c Reminder: this comes from an antisymmetrix matrix 2380c in case we want to pull back the 2d munu-matrix 2381c g_munu_rot, rotated AO matrix for dia part 2382c storing only diag + off-diag elements 2383c Reminder: this comes from a symmetric matrix 2384c in case we want to pull back the 2d munu-matrix 2385 call reorder_munu(g_munu_rot ,nat,nlst,nbf,Ndir_munu) ! reoder-munu matrix 2386 call reorder_munu(g_munu_rot1,nat,nlst,nbf,Ndir_munu) ! reoder-munu matrix 2387 call reorder_munu(g_munu_rot2,nat,nlst,nbf,Ndir_munu) ! reoder-munu matrix 2388c ------ destroy unnecessary GAs 2389 if (.not. ga_destroy(g_munuCSdia)) call errquit( 2390 & 'create_munu4nbo-1: ga_destroy failed ',0, GA_ERR) 2391 if (.not. ga_destroy(g_munuCSpar1)) call errquit( 2392 & 'create_munu4nbo-1: ga_destroy failed ',0, GA_ERR) 2393 if (.not. ga_destroy(g_munuCSHpar)) call errquit( 2394 & 'create_munu4nbo-2: ga_destroy failed ',0, GA_ERR) 2395 if (.not. ga_destroy(g_munuCSHpar2d)) call errquit( 2396 & 'create_munu4nbo-7a: ga_destroy failed ',0, GA_ERR) 2397 if (.not. ga_destroy(g_munuCSpar12d)) call errquit( 2398 & 'create_munu4nbo-7: ga_destroy failed ',0, GA_ERR) 2399 if (.not. ga_destroy(g_tnp)) call errquit( 2400 & 'create_munu4nbo-5: ga_destroy failed ',0, GA_ERR) 2401 if (.not. ga_destroy(g_acc)) call errquit( 2402 & 'create_munu4nbo-6: ga_destroy failed ',0, GA_ERR) 2403 if (.not. ga_destroy(g_acc3)) call errquit( 2404 & 'create_munu4nbo-6: ga_destroy failed ',0, GA_ERR) 2405 if (.not. ga_destroy(g_tnp1)) call errquit( 2406 & 'create_munu4nbo-5a: ga_destroy failed ',0, GA_ERR) 2407 if (.not. ga_destroy(g_acc1)) call errquit( 2408 & 'create_munu4nbo-6a: ga_destroy failed ',0, GA_ERR) 2409 if (.not. ga_destroy(g_dens)) call errquit( 2410 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 2411 if (.not. ga_destroy(g_dens1)) call errquit( 2412 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 2413 if (.not. ga_destroy(g_c1)) call errquit( 2414 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 2415 do i=1,npol 2416 911 if (.not.ga_destroy(vectors(i))) call 2417 & errquit('create_munu4nbo: ga_destroy failed vectors', 2418 & 0,GA_ERR) 2419 enddo 2420 call wshldfile(rtdb, 2421 & g_munu_rot, ! dia term 2422 & g_munu_rot2, ! 1st term in para g_munuEPRpar1 2423 & g_acc2, ! perturbed MO vector x,y,z 2424 & g_munu_rot1, ! perturbed AO operator g_munuEPRHpar x,y,z 2425 & nlst,npol_munu, 2426 & Ndir_munu, ! OUTPUT: used in wgshiftfile(rtdb) 2427 & Natoms_munu, 2428 & atmnr_munu) 2429 if (.not. ga_destroy(g_munu_rot)) call errquit( 2430 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 2431 if (.not. ga_destroy(g_munu_rot1)) call errquit( 2432 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 2433 if (.not. ga_destroy(g_munu_rot2)) call errquit( 2434 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 2435 if (.not. ga_destroy(g_acc2)) call errquit( 2436 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 2437 return 2438 end 2439c -------------- for shieldings NMLO analysis ----------------- END 2440