1 subroutine hnd_hyperfine_zora(rtdb,basis,geom) 2c $Id$ 3 implicit none 4#include "errquit.fh" 5#include "global.fh" 6#include "mafdecls.fh" 7#include "msgids.fh" 8#include "geom.fh" 9#include "rtdb.fh" 10#include "bas.fh" 11#include "stdio.fh" 12#include "apiP.fh" 13#include "prop.fh" 14#include "bgj.fh" 15#include "case.fh" 16#include "zora.fh" 17 integer rtdb ! [input] rtdb handle 18 integer basis ! [input] basis handle 19 integer geom ! [input] geometry handle 20 integer nclosed(2), nopen(2), nvirt(2), ndens, nbf, nmo 21 integer nat,nat_slc,typeprop, 22 & ixy,ix,iy,iz,iatom,iocc,ifld,ioff 23 integer alo(3),ahi(3),blo(3),bhi(3),clo(3),chi(3) 24 integer dlo(3), dhi(3) 25 integer l_occ,k_occ,l_eval,k_eval 26 integer l_dia,k_dia,l_para,k_para,l_tot,k_tot 27 integer l_xyz,k_xyz,l_zan,k_zan 28 integer l_cfine,k_cfine 29 integer l_AtNr,k_AtNr 30 integer icalczora,type_nmrdata 31 integer g_dens(3),type_NMR, 32 & g_d1, g_rhs, g_rhs0, g_u 33 integer ga_dia_hfine,ga_para1, 34 & ga_h01_hfine,ga_Fji_hfine,g_AtNr1 35 integer vectors(2), geomnew, i, j, ij, g_xc(3) 36 double precision atn, tol2e, val 37 double precision a(6),axs(3,3),eig(3),conv2MHz 38 logical dft_zoraNMR_read 39 character*255 zorafilename 40 character*3 scftyp 41 character*16 tag 42 character*32 element 43 character*256 cphf_rhs, cphf_sol 44 character*2 symbol 45 integer ld(2),cbuf,ind 46 logical cphf2, file_write_ga, file_read_ga, cphf 47 external cphf2, file_write_ga, file_read_ga, cphf 48 external giao_aotomo,mat_transpose ! FA-09-16-10 49 double precision val1,coeffpol,ac_occ(2),small,trace, 50 & isotr,aniso, 51 & span,skew,asym,rtemp 52 data small /1.0d-10/ 53 logical oskel, status 54 data tol2e /1.0d-10/ 55 integer ic,npol 56c ----- Definitions for NLMO analysis ---- START 57 integer i1,j1,acc_vec,l_tvec,k_tvec,hypfile 58 integer g_tvec ! for nbo analysis 59c ----- Definitions for NLMO analysis ---- END 60 external get_par_hfine,get_dia_hfine, 61 & add_H10_hfine,skip_cphf,get_P10, 62 & get_convfactor_hfine,sort_asc, 63 & dft_zoraNMR_read,create_munu4nbo_hyp, 64 & get_densZ4, 65 & dft_zoraCPHF_write,dft_zoraCPHF_read 66 67 integer ntot,ispin,indA,indB,nocc(2),nind_jk 68 integer debug_giaox 69 logical skip_cphf_ev_hyp 70 71 debug_giaox=0 ! =1 for debugging print outs of matrices 72 if (ga_nodeid().eq.0) then 73 write(LuOut,*) 74 call util_print_centered(LuOut, 75 $ 'Hyperfine Tensor (in au)', 23, .true.) 76 write(LuOut,*) 77 end if 78c 79c Current CPHF does not handle symmetry 80c Making C1 geometry and store it on rtdb 81c 82c Integral initialization 83 call int_init(rtdb,1,basis) 84 call schwarz_init(geom,basis) 85 call hnd_giao_init(basis,1) 86 call scf_get_fock_param(rtdb,tol2e) 87c 88c Find out from rtdb which atoms we need to calculate hyperfine for 89c Get number of atoms (all or number from rtdb) 90c Get which atoms (all or some read from rtdb) 91c Allocate arrays which will hold atomic information (k_zan and k_xyz) 92 93 status = rtdb_parallel(.true.) 94c ------- Read (nat,atmnr) --------- START 95 status=geom_ncent(geom,nat) 96 if (.not.ma_alloc_get( 97 & mt_int,nat,'nmt tmp',l_AtNr,k_AtNr)) 98 & call errquit('hnd_hyperfine_zora: ma failed',0,MA_ERR) 99 typeprop=3 ! =1 EFG =2 Shieldings =3 Hyperfine 100 call get_slctd_atoms(nat_slc, ! out: selected atoms 101 & int_mb(k_AtNr),! out: list of selected atom nr. 102 & nat, ! in : total nr atoms in molecule 103 & rtdb, ! in : rdt handle 104 & typeprop) ! in : =1,2,3=EFG,Shieldings,Hyperfine 105 if (ga_nodeid().eq.0) then 106 write(*,8) nat_slc 107 8 format('nat_slc=',i4) 108 do i=1,nat_slc 109 write(*,7) i,int_mb(k_AtNr+i-1) 110 7 format(' In hnd_hyperfine_zora:: atomnr(',i3,')=',i5) 111 enddo 112 endif 113c ------- Read (nat,atmnr) --------- END 114 if (.not. ma_push_get(mt_dbl,nat_slc,'nmr hfine',l_cfine,k_cfine)) 115 & call errquit('hnd_hfine: ma_push_get failed k_cfine',0,MA_ERR) 116 if (.not. ma_push_get(mt_dbl,3*nat_slc,'nmr at',l_xyz,k_xyz)) 117 & call errquit('hnd_hfine: ma_push_get failed k_xyz',0,MA_ERR) 118 if (.not. ma_push_get(mt_dbl,nat_slc,'nmr zan',l_zan,k_zan)) 119 & call errquit('hnd_hfine: ma_push_get failed k_zan',0,MA_ERR) 120 hypfile=0 ! not doing NLMO analysis by default 121 status=rtdb_get(rtdb,'prop:hypfile',mt_int,1,hypfile) ! for NLMO analysis 122 if (hypfile.eq.1) then ! ------- hypfile-if++++ START 123c Allocate memory for l_tvec,k_tvec --- start 124 if (.not.ma_alloc_get(mt_dbl,nat_slc*3*3,'tvec',l_tvec,k_tvec)) 125 & call errquit('hnd_hfine: ma failed',0,MA_ERR) 126 call dcopy(nat_slc*3*3,0.0d0,0,dbl_mb(k_tvec),1) ! reset 127c Allocate memory for l_tvec,k_tvec --- end 128 endif ! ------------------------ hypfile-if++++ END 129 130 do ixy = 0, nat_slc-1 131 if (.not. geom_cent_get(geom, int_mb(k_AtNr+ixy), tag, 132 & dbl_mb(k_xyz+3*ixy),dbl_mb(k_zan+ixy))) 133 & call errquit('hnd_hfine: geom_cent_tag failed',0, GEOM_ERR) 134 enddo 135 136 call get_convfactor_hfine(geom, ! in : geom handle 137 & dbl_mb(k_cfine),! out : conversion factor for hyperfine calc. 138 & dbl_mb(k_zan), 139 & nat_slc,int_mb(k_AtNr)) 140c Get Unperturbed MO vectors and eigenvalues 141c First allocate some memory for occupation numbers and eigenvalues 142 143 if (.not. bas_numbf(basis,nbf)) call 144 & errquit('hnd_hfine: could not get nbf',0, BASIS_ERR) 145c ++++++ Reading hyperfine data from file ++++++ START 146c Note.- lbl_nmrhyp defined in zora.fh 147 call util_file_name(lbl_nmrhyp,.false.,.false.,zorafilename) 148 icalczora = 0 ! initialize the flag 149 type_nmrdata=2 ! =1,2,3=shieldings,hyperfine,gshift 150 if (.not.dft_zoraNMR_read(zorafilename, 151 & type_nmrdata, 152 & nbf,nat_slc, 153 & g_AtNr1, 154 & ga_dia_hfine, 155 & ga_para1, 156 & ga_h01_hfine, 157 & ga_Fji_hfine)) icalczora=1 158c Note.- If I print the GAs here it gets freezed 159c ++++++ Reading hyperfine data from file ++++++ END 160 if (.not. ma_push_get(mt_dbl,2*nbf,'occ num',l_occ,k_occ)) call 161 & errquit('hnd_hfine: ma_push_get failed k_occ',0,MA_ERR) 162 if (.not. ma_push_get(mt_dbl,2*nbf,'eigenval',l_eval,k_eval)) call 163 & errquit('hnd_hfine: ma_push_get failed k_eval',0,MA_ERR) 164 call hnd_prp_vec_read(rtdb,geom,basis,nbf,nclosed,nopen, 165 & nvirt,scftyp,vectors,dbl_mb(k_occ), 166 & dbl_mb(k_eval),nmo) 167c ------ define npol ----- START 168 if (.not. rtdb_get(rtdb, 'dft:ipol', mt_int, 1, npol)) then 169 if (scftyp .eq. 'UHF') then 170 npol=2 171 else if (scftyp .eq. 'RHF') then 172 npol=1 173 endif 174 endif 175c ------ define npol ----- END 176c if (ga_nodeid().eq.0) 177c & write(*,*) 'npol=',npol 178 179 if (npol.eq.1) then 180 nocc(1)=nclosed(1) 181 nocc(2)=0 182 else if (npol.eq.2) then 183 nocc(1)=nopen(1) 184 nocc(2)=nopen(2) 185 endif 186 if (debug_giaox.eq.1) then 187 write(*,10) nocc(1),nocc(2), 188 & nopen(1),nopen(2), 189 & nclosed(1),nclosed(2), 190 & nvirt(1),nvirt(2),scftyp 191 10 format('nocc =(',i3,',',i3,') ', 192 & 'nopen=(',i3,',',i3,') ', 193 & 'nclos=(',i3,',',i3,') ', 194 & 'nvirt=(',i3,',',i3,') scftyp=',a,')') 195 endif 196c ---- calculate total occupation alpha and beta --- START 197 do i=1,npol 198 val=0.0d0 199 ind=nbf*(i-1) 200 do j=1,nocc(i) 201 val=val+dbl_mb(k_occ+ind+j-1) 202 enddo 203 ac_occ(i)=val 204 enddo 205c if (ga_nodeid().eq.0) then 206c write(*,2) ac_occ(1),ac_occ(2),nocc(1),nocc(2) 207c 2 format('In hnd_hfine: ac_occ=(', 208c & f15.8,',',f15.8,') nocc=(', 209c & i4,',',i4,')') 210c endif 211c ---- calculate total occupation alpha and beta --- END 212c ----- calc coeff_pol --------- START 213 if (ac_occ(1) .ne. ac_occ(2)) then 214 coeffpol=1.0d0/(ac_occ(1)-ac_occ(2)) 215 else 216 write(*,1) nocc(1),nocc(2) 217 1 format('Error in hnd_hyperfine_zora(): ', 218 & 'noc=(',i3,',',i3,') ', 219 & '-> closed shell system not allowed!') 220 stop 221 endif 222c if (ga_nodeid().eq.0) 223c & write(*,*) 'coeffpol=',coeffpol 224c ----- calc coeff_pol --------- END 225 226c ------ Store nopen in rtdb so that CPHF routine is happy ---- START 227c WARNING: For restricted calc nocc(1)=9 nocc(2)=0 <=== PROBLEM 228 if (npol.eq.2) then 229 if (.not. rtdb_put(rtdb, 'scf:nopen', 230 & MT_INT, 1, nocc(1)-nocc(2))) 231 * call errquit('scfinit:rtdbput nopen failed', 232 & nocc(1)-nocc(2), 233 & RTDB_ERR) 234 endif 235c ------ Store nopen in rtdb so that CPHF routine is happy ---- END 236c 237c Get Unperturbed Density Matrix 238 call hnd_prp_get_dens(rtdb,geom,basis, 239 & g_dens, ! out 240 & ndens,scftyp, 241 & nclosed,nopen,nvirt) 242 ntot=0 243 do ispin=1,npol 244 ntot=ntot+nocc(ispin)*nvirt(ispin) 245 enddo 246 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs)) 247 & call errquit('hnd_gshift: ga_create failed g_rhs',0,GA_ERR) 248 call ga_zero(g_rhs) 249 250 if (debug_giaox.eq.1) then 251 write(*,*) 'after-creating g_rhs ...' 252 write(*,102)npol,nocc(1),nocc(2), 253 & nclosed(1),nclosed(2), 254 & nvirt(1),nvirt(2),scftyp,ndens 255 102 format('BEF pre-fock::npol=',i3,' nocc=(',i3,',',i3,') ', 256 & 'nclos=(',i3,',',i3,') ', 257 & 'nvirt=(',i3,',',i3,') scftyp=',a,',ndens=',i3,')') 258 endif 259 call add_H10_hfine(g_rhs, ! out: ga_rhs(a,i)=ga_rhs(a,i)+H10(a,i) 260 & ga_Fji_hfine,vectors, 261 & basis,nbf,nmo,npol,nocc,nvirt) 262 if (debug_giaox.eq.1) then 263 if (ga_nodeid().eq.0) then 264 write(*,*) 'Aft. add_H10()' 265 write(*,*) '---- g_rhs-aft-add_H10 -------- START' 266 endif 267 call ga_print(g_rhs) 268 if (ga_nodeid().eq.0) 269 & write(*,*) '---- g_rhs-aft-add_H10 -------- END' 270 endif 271 272c ----> SKIPPING CPHF ------------------ START 273c if (ga_nodeid().eq.0) 274c & write(*,*) 'WARNING : Using skip_cphf ...' 275c call skip_cphf(g_rhs, ! IN/OUT 276c & dbl_mb(k_eval), ! IN: energy eigenvalues 277c & nocc,nvirt,nmo,npol) 278c ----> SKIPPING CPHF ------------------ END 279 280c -------free allocated memory -------------- START 281 if (.not.ma_pop_stack(l_eval)) call 282 & errquit('hnd_gshift: ma_pop_stack failed k_eval',0,MA_ERR) 283 if (.not.ma_pop_stack(l_occ)) call 284 & errquit('hnd_gshift: ma_pop_stack failed k_occ',0,MA_ERR) 285c -------free allocated memory -------------- END 286 if (debug_giaox.eq.1) then 287 if (ga_nodeid().eq.0) 288 & write(*,*) '---- Reading INPUT-cphf: g_rhs -------- START' 289 call ga_print(g_rhs) 290 if (ga_nodeid().eq.0) 291 & write(*,*) '---- Reading INPUT-cphf: g_rhs -------- END' 292 endif 293c ------- creating dummy var -------- START 294 ntot=0 295 do ispin=1,npol 296 ntot=ntot+nocc(ispin)*nocc(ispin) 297 enddo 298 if(.not.ga_create(MT_DBL,ntot,3,'RHS',-1,-1,g_rhs0)) 299 & call errquit('get_prelim_fock: ga_create failed g_rhs0', 300 & 0,GA_ERR) 301 call ga_zero(g_rhs0) 302c ------- creating dummy var -------- END 303c if (ga_nodeid().eq.0) 304c & write(*,*) 'WARNING: SKIP cphf ...' 305c goto 123 306 if(.not.rtdb_get(rtdb,'zora:skip_cphf_ev_hyp', 307 & mt_log,1,skip_cphf_ev_hyp)) 308 & skip_cphf_ev_hyp = .false. 309 if (.not.(skip_cphf_ev_hyp)) then ! ----- if-skip-cphf--START 310c if (ga_nodeid().eq.0) 311c & write(*,*) 'COMPUTE cphf hyp data ...' 312c Write ga_rhs to disk 313 call cphf_fname('cphf_rhs',cphf_rhs) 314 call cphf_fname('cphf_sol',cphf_sol) 315 if(.not.file_write_ga(cphf_rhs,g_rhs)) call errquit 316 $ ('hnd_hfine: could not write cphf_rhs',0, DISK_ERR) 317 call schwarz_tidy() 318 call int_terminate() 319c 320c Call the CPHF routine 321c 322c We do need to tell the CPHF that the density is skew symmetric. 323c Done via rtdb, put cphf:skew .false. on rtdb and later remove it. 324 if (.not. rtdb_put(rtdb, 'cphf:skew', mt_log, 1,.false.)) call 325 $ errquit('hnd_hfine: failed to write skew ', 0, RTDB_ERR) 326 if (.not.cphf2(rtdb)) call errquit 327 $ ('hnd_hfine: failure in cphf ',0, RTDB_ERR) 328 if (.not. rtdb_delete(rtdb, 'cphf:skew')) call 329 $ errquit('hnd_hfine: rtdb_delete failed ', 0, RTDB_ERR) 330c 331c Occ-virt blocks are the solution pieces of the CPHF 332c Read solution vector from disk and put solutions in U matrices 333 call ga_zero(g_rhs) 334 if(.not.file_read_ga(cphf_sol,g_rhs)) call errquit 335 $ ('hnd_hfine: could not read cphf_rhs',0, DISK_ERR) 336c ----- write CPHF data to file ----------------- START 337c Note.- lbl_cphfhyp defined in zora.fh 338 call util_file_name(lbl_cphfhyp,.false.,.false.,zorafilename) 339c if (ga_nodeid().eq.0) 340c & write(*,*) '---------- g_rhs-hyp --------- START' 341c call ga_print(g_rhs) 342c if (ga_nodeid().eq.0) 343c & write(*,*) '---------- g_rhs-hyp --------- END' 344 call dft_zoraCPHF_write( 345 & zorafilename, ! in: filename 346 & npol, ! in: nr polarization 347 & nocc, ! in: nr occupied MOs 348 & nvirt, ! in: nr virtual MOs 349 & nbf, ! in: nr basis functions 350 & vectors, ! in: MOs 351 & g_rhs0, ! out: (ntot,3) GA matrix 352 & g_rhs) ! in: (nocc*nvirt,3) GA matrix 353c ----- write CPHF data to file ----------------- END 354 else 355 call ga_zero(g_rhs) 356 call ga_zero(g_rhs0) 357 do i=1,npol 358 call ga_zero(vectors(i)) 359 enddo 360 if (ga_nodeid().eq.0) 361 & write(*,*) 'READ cphf hyperfine data from file ...' 362 call util_file_name(lbl_cphfhyp,.false.,.false.,zorafilename) 363 call dft_zoraCPHF_read( 364 & zorafilename, ! in: filename 365 & npol, ! in: nr polarization 366 & nocc, ! in: nr occupied MOs 367 & nvirt, ! in: nr virtual MOs 368 & nbf, ! in: nr basis functions 369 & vectors, ! out: MOs 370 & g_rhs0, ! out: (ntot,3) GA matrix 371 & g_rhs) ! out: (nocc*nvirt,3) GA matrix 372c if (ga_nodeid().eq.0) 373c & write(*,*) '---------- g_rhs-hyp-read --------- START' 374c call ga_print(g_rhs) 375c if (ga_nodeid().eq.0) 376c & write(*,*) '---------- g_rhs-hyp--read --------- END' 377 endif ! ----- if-skip-cphf--END 378 379 123 continue ! SKIPPING cphf 380 381 if (debug_giaox.eq.1) then 382 if (ga_nodeid().eq.0) 383 & write(*,*) '---- Reading sol-cphf: g_rhs -------- START' 384 call ga_print(g_rhs) 385 if (ga_nodeid().eq.0) 386 & write(*,*) '---- Reading sol-cphf: g_rhs -------- END' 387 endif 388 type_NMR=2 ! =1,2,3=shieldings,hyperfine,gshift 389 call get_P10(g_d1, ! out: Perturbed density matrix occ-occ contrib 390 & type_NMR, ! in: =1,2,3=shieldings,hyperfine,gshift 391 & g_rhs, ! in: accumulated rhs expression 392 & g_rhs0, ! in: from get_prelim_fock() 393 & vectors,g_CiFull, 394 & nbf,nmo,npol,nocc,nvirt, 395 & do_zora,do_NonRel,not_zora_scale,rtdb, 396 & lbl_nlmohyp) 397 if (.not.ga_destroy(g_rhs)) call 398 & errquit('hnd_hfine_zora: ga_destroy failed g_rhs',0,GA_ERR) 399 if (.not.ga_destroy(g_rhs0)) call 400 & errquit('hnd_hfine_zora: ga_destroy failed g_rhs0',0,GA_ERR) 401 402c Now we have in g_d1(nmo,nmo,3) the derivative densities and 403c hence we can calculate the contributions to the hyperfine tensor 404 if (.not. ma_push_get(mt_dbl,9*nat_slc,'sh para',l_para,k_para)) 405 & call errquit('hnd_hfine: ma_push_get failed k_para',0,MA_ERR) 406 if (.not. ma_push_get(mt_dbl,9*nat_slc,'sh dia',l_dia,k_dia)) call 407 & errquit('hnd_hfine: ma_push_get failed k_dia',0,MA_ERR) 408 if (.not. ma_push_get(mt_dbl,9*nat_slc,'sh para',l_tot,k_tot)) 409 & call errquit('hnd_hfine: ma_push_get failed k_tot',0,MA_ERR) 410c Before we start getting the integrals we need to reinitialize the 411c integrals. They were terminated by the cphf. 412 call int_init(rtdb,1,basis) 413 call schwarz_init(geom,basis) 414 call hnd_giao_init(basis,1) 415 call get_par_hfine(dbl_mb(k_para), ! OUT: para-hfine tensor 416 & ga_h01_hfine, ! IN : h01 matrix 417 & coeffpol, 418 & g_d1, ! IN : Perturbed density matrix 419 & basis,nbf,nat_slc,rtdb) 420 call get_dia_hfine(dbl_mb(k_dia), ! OUT: dia-hfine tensor 421 & ga_dia_hfine, ! IN 422 & coeffpol, 423 & basis,nbf,nat_slc,rtdb) 424c +++++++++ print-total-pardia-transferred +++ START 425 if (ga_nodeid().eq.0) 426 & write(*,105) nat_slc 427 105 format('NATOMS=',i3) 428 ic=1 429 do iatom = 1, nat_slc 430 do ix = 1, 3 431 do iy = 1, 3 432 if (ga_nodeid().eq.0) then 433 write(*,19) ix,iy,iatom, 434 & dbl_mb(k_dia +ic-1), 435 & dbl_mb(k_para+ic-1), 436 & dbl_mb(k_dia +ic-1)+ 437 & dbl_mb(k_para+ic-1) 438 19 format('NW:(dia,par,dia+par)(', 439 & i1,',',i1,',',i1,')=(', 440 & f15.8,' ',f15.8,' ',f15.8,' )') 441 endif 442 ic=ic+1 443 enddo ! end-loop-iy 444 enddo ! end-loop-ix 445 enddo ! end-loop-iatom 446c +++++++++ print-total-dia-transferred +++ END 447 do ispin=1,ndens 448 if (.not.ga_destroy(g_dens(ispin))) call 449 & errquit('hnd_hfine_zora: ga_destroy failed g_dens', 450 & 0,GA_ERR) 451 enddo 452c 453c Print out tensor information, and write to Ecce file if necessary 454 status = rtdb_parallel(.false.) 455 if (ga_nodeid().gt.0) goto 300 456 acc_vec=0 ! For NMLO analysis 457 call ecce_print_module_entry('nmr') 458 do iatom = 1, nat_slc 459 ioff = (iatom-1)*9 460 if (.not. geom_cent_get(geom, int_mb(k_AtNr+iatom-1), tag, 461 & dbl_mb(k_xyz), dbl_mb(k_zan))) call 462 & errquit('hnd_hfine: geom_cent_tag failed',0,GEOM_ERR) 463 if (.not. geom_tag_to_element(tag, symbol, element, atn)) call 464 & errquit('hnd_hfine: geom_tag_to_element failed',0,GEOM_ERR) 465c 466c Print tensor pieces and sum for total hyperfine tensor 467 if (ga_nodeid().eq.0) then 468c---- 2nd format: 2col hyperfine tensor (1scol: au 2ndcol: MHz) ---- START 469 conv2MHz=dbl_mb(k_cfine+iatom-1) 470 write(luout,9700) int_mb(k_AtNr+iatom-1),symbol ! Showing original atm nr. 471 write(luout,97) 472 write(luout,112) ' FC, SD terms:' 473 trace=(dbl_mb(k_dia+ioff ) + 474 & dbl_mb(k_dia+ioff+4) + 475 & dbl_mb(k_dia+ioff+8))/3.0d0 476 write(luout,99) trace,trace*conv2MHz 477 99 format(' isotropic',t16,'A(au)', 478 & f12.4,t60,'A(MHz)',f12.4) 479 do iy=1,3 480 iz=3*(iy-1) 481 write (luout,98) (dbl_mb(k_dia+ioff+ix-1),ix=iz+1,iz+3), 482 & (dbl_mb(k_dia+ioff+ix-1)*conv2MHz, 483 & ix=iz+1,iz+3) 484 enddo 485 write(luout,112) ' PSO-Spin-Orbit terms:' 486 trace=(dbl_mb(k_para+ioff ) + 487 & dbl_mb(k_para+ioff+4) + 488 & dbl_mb(k_para+ioff+8))/3.0d0 489 write(luout,99) trace,trace*conv2MHz 490 do iy=1,3 491 iz=3*(iy-1) 492 write (luout,98) (dbl_mb(k_para+ioff+ix-1),ix=iz+1,iz+3), 493 & (dbl_mb(k_para+ioff+ix-1)*conv2MHz, 494 & ix=iz+1,iz+3) 495 enddo 496 do ix = 0, 8 497 dbl_mb(k_tot+ioff+ix) = dbl_mb(k_dia+ioff+ix) + 498 & dbl_mb(k_para+ioff+ix) 499 enddo 500 write(luout,112) ' Total hyperfine coupling tensor:' 501 trace=(dbl_mb(k_tot+ioff ) + 502 & dbl_mb(k_tot+ioff+4) + 503 & dbl_mb(k_tot+ioff+8))/3.0d0 504 write(luout,99) trace,trace*conv2MHz 505 do iy=1,3 506 iz=3*(iy-1) 507 write (luout,98) (dbl_mb(k_tot+ioff+ix-1),ix=iz+1,iz+3), 508 & (dbl_mb(k_tot+ioff+ix-1)*conv2MHz, 509 & ix=iz+1,iz+3) 510 enddo 511c---- 2nd format: 2col hyperfine tensor (1scol: au 2ndcol: MHz) ---- END 512 endif 513 514 do ix = 0, 8 515 dbl_mb(k_para+ioff+ix) = dbl_mb(k_dia+ioff+ix) + 516 & dbl_mb(k_para+ioff+ix) 517 enddo 518c 519c Print total hyperfine tensor 520 521 if (ga_nodeid().eq.0) then 522c write(luout,9802) (dbl_mb(k_para+ioff+ix-1),ix=1,9) 523c 524c Diagonalize total tensor 525c Order in a: xx xy yy xz yz zz 526 a(1) = dbl_mb(k_para+ioff) 527 a(2) = dbl_mb(k_para+ioff+1) 528 a(3) = dbl_mb(k_para+ioff+4) 529 a(4) = dbl_mb(k_para+ioff+2) 530 a(5) = dbl_mb(k_para+ioff+5) 531 a(6) = dbl_mb(k_para+ioff+8) 532 ij = 0 533 do 241 i = 1, 3 534 do 241 j = 1, i 535 ij = ij + 1 536 axs(i,j) = a(ij) 537 axs(j,i) = a(ij) 538 241 continue 539 call hnd_diag(axs,eig,3,.true.,.true.) 540 call sort_asc(eig,axs) ! sort in ascending order 541 isotr =(eig(1) + eig(2) + eig(3))/3.0d0 542 aniso = eig(1) -(eig(2) + eig(3))/2.0d0 543c ++++++++++ get eigenvectors for NMLO analysis +++++ START 544 hypfile=0 ! not doing NLMO analysis by default 545 status=rtdb_get(rtdb,'prop:hypfile',mt_int,1,hypfile) ! for NLMO analysis 546 if (hypfile.eq.1) then 547 do i1=1,3 548 do j1=1,3 549 dbl_mb(k_tvec+acc_vec)=axs(i1,j1) 550 acc_vec=acc_vec+1 551 enddo 552 enddo 553 endif 554c ++++++++++ get eigenvectors for NMLO analysis +++++ END 555c ---- Calc: span,skew,asymm ----- START 556 span = eig(3) - eig(1) 557 skew = 0.0d0 558 if (abs(span)>small) 559 & skew = 3.0d0*(isotr-eig(2))/span 560 rtemp = eig(3) - isotr 561 asym = 0.0d0 562 if (abs(rtemp)>small) 563 & asym = (eig(2)-eig(1))/rtemp 564 write(luout,101) 565 101 format(/1x,'In principal axis representation:', 566 & ' total span,skew and asymm') 567 write(luout,100) span,skew,asym, 568 & span*conv2MHz, 569 & skew*conv2MHz, 570 & asym*conv2MHz 571 100 format(t9,'span',t21,'skew',t32,'asymm',t52, 572 & 'span',t64,'skew',t75,'asymm'/3f12.4, 573 & t44,3f12.4/) 574c ---- Calc: span,skew,asymm ----- END 575 write(luout,9986) (ix,ix=1,3),(ix,ix=1,3) 576 write(luout,4489) (eig(ix),ix=1,3), 577 & (eig(ix)*conv2MHz,ix=1,3) 578 4489 format(3f12.4,t44,3f12.4,/) 579 do iy=1,3 580 write(luout,9983) iy,axs(iy,1),axs(iy,2),axs(iy,3) 581 enddo 582 write(luout,'(//)') 583c 584c Print Ecce information 585 586 call ecce_print1_char('atom name',symbol,1) 587 call ecce_print2('hyperfine tensor',MT_DBL, 588 & dbl_mb(k_para+ioff),3,3,3) 589 call ecce_print1('hyperfine isotropic',MT_DBL,isotr,1) 590 call ecce_print1('hyperfine anisotropy',MT_DBL,aniso,1) 591 call ecce_print1('hyperfine eigenvalues',MT_DBL,eig,3) 592 call ecce_print2('hyperfine eigenvectors',MT_DBL,axs, 593 & 3,3,3) 594 endif 595 enddo 596 call ecce_print_module_exit('nmr','ok') 597300 call ga_sync() 598 status = rtdb_parallel(.true.) 599 600c ---- Destroy stored ga arrays ------ START 601 if (.not. ga_destroy(g_d1)) call errquit( 602 & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 603 if (.not. ga_destroy(ga_dia_hfine)) call errquit( 604 & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 605 if (.not. ga_destroy(ga_h01_hfine)) call errquit( 606 & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 607 if (.not. ga_destroy(ga_Fji_hfine)) call errquit( 608 & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 609c ---- Destroy stored ga arrays ------ END 610 611 hypfile=0 ! not doing NLMO analysis by default 612 status=rtdb_get(rtdb,'prop:hypfile',mt_int,1,hypfile) ! for NLMO analysis 613 if (hypfile.eq.1) then ! ------- hypfile-if++++ START 614 if (.not. ga_create(mt_dbl,1,nat_slc*9, 615 & 'munu4nbo: g_tvec',0,0,g_tvec)) 616 $ call errquit('munu4nbo: g_tvec', 0,GA_ERR) 617 call ga_dgop(msg_efgs_col,dbl_mb(k_tvec),nat_slc*9,'+') 618 call ga_put(g_tvec,1,1,1,nat_slc*9,dbl_mb(k_tvec),1) 619 call create_munu4nbo_hyp(rtdb,g_tvec, 620 & coeffpol, 621 & nat_slc,int_mb(k_AtNr), 622 & basis,npol,nocc,nvirt,nmo, 623 & dbl_mb(k_cfine)) 624 if (.not. ga_destroy(g_tvec)) call errquit( ! destroy GA 625 & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 626 if (.not.ma_free_heap(l_tvec)) call 627 & errquit('hnd_hyperfine_zora: ma_free_heap l_tvec',0,MA_ERR) 628 endif ! ------------------------ hypfile-if++++ END 629c ---- Destroy stored ga arrays ------ START 630c if (.not. ga_destroy(g_d1)) call errquit( 631c & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 632c if (.not. ga_destroy(ga_dia_hfine)) call errquit( 633c & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 634c if (.not. ga_destroy(ga_h01_hfine)) call errquit( 635c & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 636c if (.not. ga_destroy(ga_Fji_hfine)) call errquit( 637c & 'hnd_hyperfine_zora: ga_destroy failed ',0, GA_ERR) 638c ---- Destroy stored ga arrays ------ END 639c 640c Clean up all remaining memory 641 if (.not.ma_pop_stack(l_tot)) call 642 & errquit('hnd_hyperfine_zora: ma_pop_stack failed k_tot', 643 & 0,MA_ERR) 644 if (.not.ma_pop_stack(l_dia)) call 645 & errquit('hnd_hyperfine_zora: ma_pop_stack failed k_dia', 646 & 0,MA_ERR) 647 if (.not.ma_pop_stack(l_para)) call 648 & errquit('hnd_hyperfine_zora: ma_pop_stack failed k_para', 649 & 0,MA_ERR) 650 do i=1,npol 651 911 if (.not.ga_destroy(vectors(i))) call 652 & errquit('hnd_hyperfine_zora: ga_destroy failed vectors', 653 & 0,GA_ERR) 654 enddo 655 if (.not. ga_destroy(g_AtNr1)) call errquit( 656 & 'hnd_giaox_zora: ga_destroy failed ',0, GA_ERR) 657 if (.not.ma_pop_stack(l_zan)) call 658 & errquit('hnd_hyperfine_zora: ma_pop_stack failed k_zan', 659 & 0,MA_ERR) 660 if (.not.ma_pop_stack(l_xyz)) call 661 & errquit('hnd_hyperfine_zora: ma_pop_stack failed k_xyz', 662 & 0,MA_ERR) 663 if (.not.ma_free_heap(l_AtNr)) call 664 & errquit('hnd_hyperfine_zora: ma_free_heap l_tmp', 665 & 0, MA_ERR) 666 if (.not.ma_pop_stack(l_cfine)) call 667 & errquit('hnd_hyperfine_zora: ma_pop_stack failed l_cfine', 668 & 0,MA_ERR) 669 call schwarz_tidy() 670 call int_terminate() 671 return 672 673 97 format(/' ----- A-tensor in atomic units ----',t44, 674 1 ' ------------- A (MHz) ------------') 675 98 format(3f12.4,t44,3f12.4) ! for showing 2col hyperfine tensor 676 112 format(/a) 677 7000 format(/,10x,'NMR hyperfine cannot be calculated for UHF', 678 1 ' or ROHF wave functions: use RHF') 679 9700 format(6x,'Atom: ',i4,2x,a2) 680 9800 format(8x,'Diamagnetic',/,3(3F12.4,/)) 681 9801 format(8x,'Paramagnetic',/,3(3F12.4,/)) 682 9802 format(8x,'Total hyperfine Tensor',/,3(3F12.4,/)) 683 9983 format(i1,1x,f10.4,f12.4,f12.4) 684c 9983 format(6x,i1,3x,3f12.4) 685 9985 format(10x,3f12.4,/) 686c 9986 format(10x,'Principal Components and Axis System',/,10x, 687c 1 3(7x,i1,4x)) 688 9986 format(1x,'Principal Components and Axis System',/,1x, 689 1 3(7x,i1,4x),t44,3(7x,i1,4x)) 690 9987 format(10x,' isotropic = ',f12.4,/, 691 1 10x,'anisotropy = ',f12.4,/) 692 9999 format( 693 1 /,10x,41(1h-),/, 694 2 10x,'Hyperfine Tensors (in au)',/, 695 3 10x,41(1h-),/) 696 end 697c 698c --------------- For NMLO analysis ------------------ START 699 subroutine get_tvec(dia, ! in : dia tensor 700 & tvec, ! out: eigenvectors 701 & acc_vec, ! in/out: counter for eigenvectors 702 & iat, ! in : atom nr. 703 & hypfile) ! in : NMLO-hyp flag 704 implicit none 705#include "errquit.fh" 706#include "global.fh" 707#include "mafdecls.fh" 708#include "msgids.fh" 709#include "rtdb.fh" 710#include "apiP.fh" 711#include "prop.fh" 712#include "bgj.fh" 713#include "bas.fh" 714#include "stdio.fh" 715 716 integer iat,acc_vec 717 double precision dia(*),tvec(*) 718 integer i,j,ic,hypfile 719 double precision vec(3,3), eig(3) 720 external hnd_diag 721 ic=9*(iat-1)+1 722 do i = 1, 3 723 do j = 1, 3 724 vec(i,j)= dia(ic) 725 ic=ic+1 726 enddo ! end-loop-i 727 enddo ! end-loop-j 728 call hnd_diag(vec,eig,3,.true.,.false.) 729c ------- copy eigenvectors, vec for NLMO analysis ----- START 730 if (hypfile.eq.1) then 731 do i=1,3 732 do j=1,3 733 tvec(acc_vec)=vec(i,j) 734 acc_vec=acc_vec+1 735 enddo 736 enddo 737 endif 738 return 739 end 740c --------------- For NMLO analysis ------------------ END 741 subroutine get_par_hfine( 742 & par, ! OUT : paramagnetic tensor 743 & ga_h01_hfine,! IN : 744 * coeffpol, ! IN : scaling factor 1/(nA-nB) 745 & g_d1, ! IN : Perturbed density matrix 746 & basis, ! IN : basis handle 747 & nbf, ! IN : nr. basis functions 748 & natoms, ! IN : nr. atoms 749 & rtdb) 750c Purpose : Assemble NMR Hyperfine: paramagnetic tensor 751c Author : Fredy Aquino 752c Date : 03-09-11 753 implicit none 754#include "errquit.fh" 755#include "global.fh" 756#include "mafdecls.fh" 757#include "msgids.fh" 758#include "rtdb.fh" 759#include "apiP.fh" 760#include "prop.fh" 761#include "bgj.fh" 762#include "bas.fh" 763#include "stdio.fh" 764#include "zora.fh" 765c Global arrays variables defined in zora.fh 766c ga_dia_hfine,ga_h01_hfine,ga_Fji_hfine 767 integer g_d1 ! input: Perturbed density matrices 768 integer ga_h01_hfine 769 integer rtdb 770 integer basis 771 integer alo(3), ahi(3), blo(3), bhi(3) 772 integer ld(2),cbuf,cbuf1,cbuf2 773 integer nbf,natoms 774 integer iatom,ix,iy,ixy,ind 775 integer debug_get_par 776 double precision val,par(*),coeffpol 777 logical oskel 778 debug_get_par=0 779 do ixy = 1, 9*natoms 780 par(ixy) = 0.0d0 ! initialize 781 enddo 782 alo(1) = 1 783 ahi(1) = nbf 784 alo(2) = 1 785 ahi(2) = nbf 786 blo(1) = 1 787 bhi(1) = nbf 788 blo(2) = 1 789 bhi(2) = nbf 790 ixy = 0 791 blo(3) = 0 792 bhi(3) = 0 793 do iatom = 1, natoms 794 do iy = 1, 3 795 blo(3) = blo(3) + 1 796 bhi(3) = bhi(3) + 1 797 do ix = 1, 3 798 alo(3) = ix 799 ahi(3) = ix 800 ixy = ixy + 1 801 val = nga_ddot_patch(g_d1 ,'n',alo,ahi, 802 & ga_h01_hfine,'n',blo,bhi)* 803 * coeffpol* 804 & (-2.0d0) ! factor from (+2) g_N B_N/(nA-nB) 805 ! (-1) from h_01 operator 806 cbuf=9*(iatom-1)+3*(ix-1)+iy ! transpose 807 par(cbuf)=val 808 enddo ! end-loop-ix 809 enddo ! end-loop-iy 810 enddo ! end-loop-iatom 811c --------- symmmetrize para ------- START 812 do iatom = 1, natoms 813 do iy = 1,3 814 do ix = iy+1,3 815 cbuf1=9*(iatom-1)+3*(ix-1)+iy ! transpose 816 cbuf2=9*(iatom-1)+3*(iy-1)+ix 817 val=(par(cbuf1)+par(cbuf2))/2.0d0 818 par(cbuf1)=val 819 par(cbuf2)=val 820 enddo ! end-loop-ix 821 enddo ! end-loop-iy 822 enddo ! end-loop-iatom 823c --------- symmmetrize para ------- END 824 return 825 end 826 827 subroutine get_dia_hfine(dia, ! OUT : dia-hfine tensor 828 & ga_dia_hfine,! IN : 829 & coeffpol, ! IN : scaling factor 830 & basis, ! IN : basis handle 831 & nbf, ! IN : nr. basis functions 832 & natoms, ! IN : nr. atoms 833 & rtdb) 834 835c Purpose : Assemble NMR Hyperfine: diamagnetic tensor 836c Author : Fredy Aquino 837c Date : 03-09-11 838 implicit none 839#include "errquit.fh" 840#include "global.fh" 841#include "mafdecls.fh" 842#include "msgids.fh" 843#include "rtdb.fh" 844#include "bas.fh" 845#include "apiP.fh" 846#include "prop.fh" 847#include "bgj.fh" 848#include "stdio.fh" 849#include "zora.fh" 850c Global arrays variables defined in zora.fh 851c --> ga_dia_hfine,ga_h01_hfine,ga_Fji_hfine 852 integer rtdb 853 integer basis ! [input] basis handle 854 integer alo(3), ahi(3), blo(3), bhi(3) 855 integer ld(2),cbuf,cbuf1,cbuf2 856 integer nbf,natoms,ispin 857 integer iatom,ix,iy,ixy 858 integer ga_dia_hfine 859 double precision dia(*),coeffpol 860 logical oskel 861 double precision val 862 do ixy = 1, 9*natoms 863 dia(ixy) = 0.0d0 ! initialize 864 enddo 865c ---- STORE: ga_dia --> dia 866 alo(1)=1 867 ahi(1)=3 868 alo(2)=1 869 ahi(2)=3 870 alo(3)=1 871 ahi(3)=natoms 872 ld(1)=3 873 ld(2)=3 874 call ga_scale(ga_dia_hfine,coeffpol) 875 call nga_get(ga_dia_hfine,alo,ahi,dia(1),ld) 876c --------- symmmetrize dia ------- START 877 do iatom = 1, natoms 878 do iy = 1,3 879 do ix = iy + 1,3 880 cbuf1=9*(iatom-1)+3*(ix-1)+iy ! transpose 881 cbuf2=9*(iatom-1)+3*(iy-1)+ix 882 val=(dia(cbuf1)+dia(cbuf2))/2.0d0 883 dia(cbuf1)=val 884 dia(cbuf2)=val 885 enddo ! end-loop-ix 886 enddo ! end-loop-iy 887 enddo ! end-loop-iatom 888c --------- symmmetrize dia ------- END 889 return 890 end 891 892 subroutine add_H10_hfine( 893 & g_rhs, ! out: accumulated rhs expression 894 & ga_Fji, ! in: Fock 1st-deriv from operator: p K x p 895 & vectors, ! in: MO coeffs 896 & basis, ! in: basis handle 897 & nbf, ! in: nr. basis functions 898 & nmo, ! in: nr. MOs (occ+virt) 899 & npol, ! in: nr. of polarizations 900 & nocc, ! in: nr. occ MOs 901 & nvirt) ! in: nr. virtual MOs 902c ============= Computing ============== 903c ga_rhs(a,i) = ga_rhs(a,i) + H10(a,i) 904c ======================================== 905 implicit none 906#include "errquit.fh" 907#include "global.fh" 908#include "mafdecls.fh" 909#include "msgids.fh" 910#include "geom.fh" 911#include "rtdb.fh" 912#include "bas.fh" 913#include "stdio.fh" 914 integer npol ! nr. of polarizations 915 logical do_zora ! logical to check if NOT doing zora calc 916 integer nbf,nmo,ispin, 917 & nocc(npol),nvirt(npol) 918 integer shift,disp 919 integer alo(3), ahi(3), 920 & blo(3), bhi(3) 921 integer vectors(npol) 922 integer g_rhs,ga_Fji 923 integer g_s10 ! scratch ga arrays 924 integer basis ! basis handle 925 logical oskel 926 external giao_aotomo 927 oskel = .false. 928c ----------- Create scratch ga-arrays ------- START 929 alo(1) = nbf 930 alo(2) = -1 931 alo(3) = -1 932 ahi(1) = nbf 933 ahi(2) = nbf 934 ahi(3) = 3*npol 935 if (.not.nga_create(MT_DBL,3,ahi,'s10 matrix', 936 & alo,g_s10)) 937 & call errquit('add_H10: nga_create failed g_s10', 938 & 0,GA_ERR) 939 call ga_zero(g_s10) 940c ----------- Create scratch ga-arrays ------- END 941c Transform H10 to MO and add to g_rhs 942c Copy: g_Fji --> g_s10 943 blo(1) = 1 944 bhi(1) = nmo 945 blo(2) = 1 946 bhi(2) = nmo 947 blo(3) = 1 948 bhi(3) = 3 949 do ispin=1,npol 950 shift=3*(ispin-1) 951 alo(1) = 1 952 ahi(1) = nmo 953 alo(2) = 1 954 ahi(2) = nmo 955 alo(3) = shift+1 956 ahi(3) = shift+3 957 call nga_copy_patch('n',ga_Fji,blo,bhi, 958 & g_s10,alo,ahi) 959 enddo ! end-loop-ispin 960 call giao_aotomo(g_s10,vectors,nocc,nvirt,npol,3,nbf) 961 do ispin=1,npol 962 shift=3*(ispin-1) 963 alo(1) = nocc(ispin)+1 964 ahi(1) = nmo 965 alo(2) = 1 966 ahi(2) = nocc(ispin) 967 alo(3) = shift+1 968 ahi(3) = shift+3 969c ------ definitions for g_rhs -------- START 970 disp=nocc(1)*nvirt(1)*(ispin-1) 971 blo(1) = disp+1 972 bhi(1) = disp+nocc(ispin)*nvirt(ispin) 973 blo(2) = 1 974 bhi(2) = 3 975c ------ definitions for g_rhs -------- END 976 call nga_add_patch(1.0d0,g_rhs,blo,bhi, 977 & 1.0d0,g_s10,alo,ahi, 978 & g_rhs,blo,bhi) 979 enddo ! end-loop-ispin 980 if (.not.ga_destroy(g_s10)) call 981 & errquit('add_H10: ga_destroy failed g_s10',0,GA_ERR) 982 return 983 end 984 subroutine get_convfactor_hfine(geom, ! in : geom handle 985 & const_hfine,! out : conversion factor for hyperfine calc. 986 & Zan, ! in : atomic numbers 987 & natoms, ! in : nr. of selected atoms 988 & atmnr) ! in : list of atom nr. indices 989 implicit none 990#include "errquit.fh" 991#include "global.fh" 992#include "mafdecls.fh" 993#include "msgids.fh" 994#include "geom.fh" 995#include "rtdb.fh" 996#include "bas.fh" 997#include "stdio.fh" 998 integer geom 999 integer iat,iat1,isonr 1000 character*16 at_tag 1001 integer natoms,atmnr(natoms) 1002 double precision gnuc,hbar,ge,emf,vl,auev,evmhz,gmhz, 1003 & pi,fac,betae,betan,conv,convf,con,gnu, 1004 & Zan(natoms) 1005 double precision const_hfine(natoms) 1006 logical status 1007 logical atom_gfac 1008 external atom_gfac 1009 data gnuc /5.05078343d-27/ ! Nuclear magneton 1010 data hbar /1.05457168d-34/ ! Planck constant over 2 pi 1011 data ge /2.002319304386d+00/ ! Electron g-factor 1012 data emf /1836.152701d+00/ ! Proton-electron mass ratio 1013 data vl /137.0359895d+00/ ! Speed of light in au 1014 data auev /27.2113961d+00/ ! Conversion from au to eV 1015 data evmhz /2.41798836d+08/ ! Conversion from eV to MHz 1016 data gmhz /2.8025d+00/ ! Conversion from Gauss to MHz 1017 1018c --- calculate constants and conversion terms ---START 1019 pi = acos(-1.0d0) 1020 fac =(4.0d0*pi/3.0d0) 1021 betae=1.0d0/(2.0d0*vl) 1022 betan=betae/emf 1023 convf=auev*evmhz 1024 con =ge*betae*betan*convf 1025c --- calculate constants and conversion terms ---END 1026c----- Allocate memory - FA 1027 do iat1 = 1, natoms 1028 iat=atmnr(iat1) 1029 if (.not. atom_gfac(Zan(iat1),gnu,isonr)) call 1030 & errquit('hnd_hfine: atom_gfac failed',0, UERR) 1031 const_hfine(iat1)=con*gnu 1032 if (ga_nodeid().eq.0) 1033 & write(*,1) iat1,iat,con,gnu,isonr,const_hfine(iat1) 1034 1 format('(con,gnu,isonr,const_hfine)(',i3,',',i3,')=(', 1035 & f15.8,',',f15.8,',',i3,',',f15.8,')') 1036 enddo ! end loop-iat1 1037 return 1038 end 1039 1040 subroutine sort_asc(eig,vec) 1041c Purpose: Sort in ascending order 1042c Taken forn hnd_diag 1043c FA-03-17-11 1044 implicit none 1045#include "errquit.fh" 1046#include "global.fh" 1047#include "msgids.fh" 1048#include "stdio.fh" 1049 integer i,j,jj 1050 double precision eig(3),vec(3,3),xx 1051 do 50 i=1,3 1052 jj=i 1053 do 30 j=i,3 1054 if( eig(j).le. eig(jj)) jj=j 1055 30 continue 1056 if(jj.eq.i) go to 50 1057 1058 xx=eig(jj) 1059 eig(jj)=eig(i) 1060 eig(i)=xx 1061 do 40 j=1,3 1062 xx=vec(j,jj) 1063 vec(j,jj)=vec(j,i) 1064 vec(j,i)=xx 1065 40 continue 1066 50 continue 1067 return 1068 end 1069c -------------------------------------------------------- 1070c -------------- for NMLO analysis ----------------- START 1071 subroutine create_munu4nbo_hyp( 1072 & rtdb, ! in: rtdb handle 1073 & g_tvec, ! in: eigenvectors or T diagonalizing matrix 1074 & coeffpol, ! in: scaling factor for hyperfine constants 1075 & nat, ! in: nr. atoms 1076 & atmnr, ! in: list of selected atoms 1077 & basis, ! in: basis handle 1078 & npol, ! in: nr. polarizaitons 1079 & nocc, ! in: nr. occ nocc(i) i=1,npol 1080 & nvirt, ! in: nr. virt nvirt(i) i=1,npol 1081 & nmo, ! in: nr. MO 1082 & const_hfine) ! in: conv2MHz 1083 1084 implicit none 1085#include "nwc_const.fh" 1086#include "errquit.fh" 1087#include "global.fh" 1088#include "bas.fh" 1089#include "mafdecls.fh" 1090#include "geom.fh" 1091#include "stdio.fh" 1092#include "rtdb.fh" 1093#include "cosmo.fh" 1094#include "msgids.fh" 1095#include "zora.fh" 1096c FA: Revised on 06-22-11 1097c ------Main outputs -------- START 1098 integer g_munu_rot, ! hyp-FCSD dia part 1099 & g_munu_rot1,g_acc2 ! hyp-PSOSO para part 1100c ------Main outputs -------- END 1101 integer rtdb,basis 1102 integer g_sdens,g_tvec 1103 integer g_munuFCSD,g_munuPSOSO, 1104 & g_tnp,g_acc, 1105 & vectors(2), 1106 & g_tnp1,g_acc1 1107 integer ispin 1108 integer nat 1109 double precision ac_val,val1,sign, 1110 & const_hfine(nat) 1111 1112 integer iat1,iat,i,j,k,m,n,ndir,ndir1 1113 integer jlo,jhi,s,nbf,nmo,nsize,nsize1,nsize2 1114 integer npol,ntot,atmnr(nat) 1115 integer ind,nlst,count,nocc(npol),nvirt(npol) 1116 integer Natoms_munu,Ndir_munu,atmnr_munu(nat) 1117c 1118c Ndir_munu, Nr. of directions stored 1119c =3 xx yy zz 1120 double precision coeff,fact,para_rot(9),coeffpol 1121 double precision tmn(2),chcdata(3) 1122 integer jlo1,jhi1,jlo2,jhi2 1123 integer g_sdens1,g_c1, 1124 & g_p10,g_munuPSOSO2d 1125 integer iind(2),jind(2),icalczora,type_NMR 1126 integer alo(3),ahi(3),elo(3),ehi(3),flo(3),fhi(3) 1127 logical dft_zoraHYP_NLMOAnalysis_read ! for read-nlmo-mat 1128 character*255 zorafilename ! for read-nlmo-mat 1129 integer arr_ind(6,2) 1130 data ((arr_ind(j,i),i=1,2),j=1,6) 1131 & /1,1,2,2,3,3,1,2,1,3,2,3/ 1132 external dft_zoraHYP_NLMOAnalysis_read,get_P10_rot, 1133 & fill_munuPSOSO_1,get_par_hfine_rot, 1134 & get2dmat,whypfile, 1135 & get_par_hfine_rot1 1136c --> To store ONLY munu principal components xx,yy,zz 1137c g_munuFCSD is created in dft_zora_Hyperfine.F 1138c size(g_munuFCSD)=nlst*ndir*nat (linear vector) 1139c g_sdens, spin density matrix 1140c nbf x nbf elements (bidimensional matrix) 1141c Legend: 1142c ndir=6 = xx, yy, zz, xy, xz, yz 1143c nbf, Nr of basis functions 1144c nlst=nbf*(nbf+1)/2 1145c nat, nr. of selected atoms out of nat 1146 1147 if (.not. bas_numbf(basis,nbf)) call errquit 1148 & ('munu: bas_numbf failed',555, BASIS_ERR) 1149 Natoms_munu=nat 1150 do i=1,Natoms_munu 1151 atmnr_munu(i)=atmnr(i) 1152 enddo 1153 Ndir_munu=3 1154 nlst=nbf*(nbf+1)/2 ! size of xx,yy,zz,xy,xz,yz chunk 1155c ++++++ Read NLMO matrices +++++++++ START 1156 ndir=6 1157 ndir1=3 1158 call util_file_name(lbl_nlmohyp,.false.,.false.,zorafilename) 1159 icalczora = 0 ! initialize the flag 1160 if (.not.dft_zoraHYP_NLMOAnalysis_read( 1161 & zorafilename, ! in : filename 1162 & nbf, ! in : nr basis functions 1163 & ndir, ! in : nr of directions: 6 = xx yy zz xy xz yz 1164 & ndir1, ! in : nr of directions: 3 = x y z 1165 & nat, ! in : list of selected atoms 1166 & nocc, ! in : nocc(i) i=1,2 nr. occupations in alpha and beta 1167 & npol, ! in: nr polarizations 1168 & g_munuFCSD, ! out: munu matrix of Fermi Contact + Spin Dipolar term 1169 & g_munuPSOSO, ! out: munu matrix of PSOSO 1170 & vectors, ! out: MOs 1171 & g_c1, ! out: perturbed MO 1172 & g_sdens)) icalczora=1 1173 call ga_scale(g_munuFCSD , coeffpol) ! scaling dia 1174 call ga_scale(g_munuPSOSO,-2.0d0*coeffpol) ! scaling para 1175 goto 10 1176 if (ga_nodeid().eq.0) 1177 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuFCSD -- START' 1178 call ga_print(g_munuFCSD) 1179 if (ga_nodeid().eq.0) 1180 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuFCSD -- END' 1181 if (ga_nodeid().eq.0) 1182 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuPSOSO -- START' 1183 call ga_print(g_munuPSOSO) 1184 if (ga_nodeid().eq.0) 1185 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_munuPSOSO -- END' 1186 if (ga_nodeid().eq.0) 1187 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_sdens ----- START' 1188 call ga_print(g_sdens) 1189 if (ga_nodeid().eq.0) 1190 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_sdens ----- END' 1191 if (ga_nodeid().eq.0) 1192 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_c1 ----- START' 1193 call ga_print(g_c1) 1194 if (ga_nodeid().eq.0) 1195 & write(*,*) 'AFT-reading-HYP-NLMO matrices: g_c1 ----- END' 1196 10 continue 1197c ++++++ Read NLMO matrices +++++++++ END 1198 call get_unique_elmat(g_sdens,g_sdens1,nlst,nbf) ! out: g_sdens1 1199 ndir=6 ! Nr. of directions: xx,yy,zz,xy,xz,yz 1200 nsize=nbf*(nbf+1)/2 ! size of xx,yy,zz,xy,xz,yz chunk 1201 nsize1=nsize*ndir ! size of whole munu per atom 1202 nsize2=nsize*ndir1 ! size of whole munu per atom 1203 if (.not. ga_create(mt_dbl,1,nsize, 1204 & 'munu4nbo: g_tnp',0,0,g_tnp)) 1205 $ call errquit('munu4nbo: g_tnp', 0,GA_ERR) 1206 call ga_zero(g_tnp) 1207 if (.not. ga_create(mt_dbl,1,nsize, 1208 & 'munu4nbo: g_tnp1',0,0,g_tnp1)) 1209 $ call errquit('munu4nbo: g_tnp1', 0,GA_ERR) 1210 call ga_zero(g_tnp1) 1211 if (.not. ga_create(mt_dbl,1,nsize, 1212 & 'munu4nbo: g_acc',0,0,g_acc)) 1213 $ call errquit('munu4nbo: g_acc', 0,GA_ERR) 1214 call ga_zero(g_acc) 1215 if (.not. ga_create(mt_dbl,1,nsize, 1216 & 'munu4nbo: g_acc1',0,0,g_acc1)) 1217 $ call errquit('munu4nbo: g_acc1', 0,GA_ERR) 1218 call ga_zero(g_acc1) 1219 alo(1) = nbf 1220 alo(2) = -1 1221 alo(3) = -1 1222 ahi(1) = nbf 1223 ntot=nocc(1)+nocc(2) 1224 ahi(2) = ntot 1225 ahi(3) = 3*nat 1226 if (.not.nga_create(MT_DBL,3,ahi,'g_acc2 matrix', 1227 & alo,g_acc2)) call 1228 & errquit('g_acc2: nga_create failed g_acc2',0,GA_ERR) 1229 call ga_zero(g_acc2) 1230 if (.not. ga_create(mt_dbl,1,nlst*3*nat, 1231 & 'munu4nbo: g_munu_rot',0,0,g_munu_rot)) 1232 $ call errquit('munu4nbo: g_munu_rot', 0,GA_ERR) 1233 if (.not. ga_create(mt_dbl,1,nlst*3*nat, 1234 & 'munu4nbo: g_munu_rot1',0,0,g_munu_rot1)) 1235 $ call errquit('munu4nbo: g_munu_rot1', 0,GA_ERR) 1236 alo(1) = nbf 1237 alo(2) = -1 1238 alo(3) = -1 1239 ahi(1) = nbf 1240 ahi(2) = nbf 1241 ahi(3) = 3 1242 if (.not.nga_create(mt_dbl,3,ahi,'g_munuPSOSO2d matrix', 1243 & alo,g_munuPSOSO2d)) call 1244 & errquit('munu4nbo: nga_create failed g_munuPSOSO2d',0,GA_ERR) 1245 if (ga_nodeid().eq.0) 1246 & write(*,*) 'CHCooooooooooooo', 1247 & ' NW-Hyperfine: Summary C+HC data [MHz] ', 1248 & 'oooooooooooooooo START' 1249 do iat1=1,nat 1250 call ga_zero(g_munuPSOSO2d) 1251 iat=atmnr(iat1) 1252 do n=1,3 ! xx,yy,zz 1253 m=n ! For principal components ONLY 1254 call ga_zero(g_acc) 1255c ----- Do: A'= T^t A T, calculate only [A']_pp --> (do n=1,3 m=n) 1256c a_pp'= \sum_i t_ip a_ii t_ip + 1257c 2 \sum_{j=2}^n \sum_{i=1}^{j-1} t_jp a_ji t_ip 1258c g_munu_rot = A' 1259c WARNING: g_munu_rot, contains several rotated matrices 1260c since the matrices are symmetric I store only 1261c the main diagonal + lower (upper) triangular 1262c matrix in a format that looks like: 1263c a_11 a_22 ... a_nn 1264c a_21 1265c a_31 a_32 1266c a_41 a_42 a_43 1267c ... 1268c a_n1 a_n2 ... a_{n(n-1)} 1269c There are two additional transformations on g_munu_rot 1270c before leaving this routine and entering wefgfile() 1271c 1. I make the diagonalized matrix traceless 1272c ===== Transform xx_munu to 2xx_munu-(yy_munu+zz_munu) = START 1273c or 3xx_munu-(xx_munu+yy_munu+zz_munu) 1274c 2. I need to do a reordering of elements so that it is 1275c compatible with wefgfile() 1276c call reorder_munu(g_munu_rot,nat,nlst,nbf,Ndir_munu) 1277c -------------------------------------------------------------- 1278 do s=1,6 1279c ------- get coeff() --- START 1280 iind(1)=1 1281 iind(2)=1 1282 jind(1)=9*(iat1-1)+3*(arr_ind(s,1)-1)+m 1283 jind(2)=9*(iat1-1)+3*(arr_ind(s,2)-1)+n 1284 call ga_gather(g_tvec,tmn,iind,jind,2) 1285 fact=1.0d0 1286 if (s.gt.3) fact=2.0d0 1287 coeff=fact*tmn(1)*tmn(2) 1288c ------- get coeff() --- END 1289c Note.- g_munuFCSD will be the (hyp-diag)_uv matrix 1290 jlo=nsize1*(iat1-1)+nsize*(s-1)+1 1291 jhi=jlo+nsize-1 1292 call ga_copy_patch('n',g_munuFCSD,1,1,jlo,jhi, 1293 & g_tnp ,1,1,1 ,nsize) 1294 call ga_add(1.0d0,g_acc,coeff,g_tnp,g_acc) 1295 enddo ! end-loop-s 1296c Note.- g_acc = \widetilde{H}_{mu nu}^{(m,m)} 1297c it is the rotated munu matrix using: T^t H T 1298 call ga_zero(g_acc1) 1299 do s=1,3 1300c ------- get coeff() --- START 1301 iind(1)=1 1302 jind(1)=9*(iat1-1)+3*(s-1)+m 1303 call ga_gather(g_tvec,tmn,iind,jind,1) 1304 1305c if (ga_nodeid().eq.0) then 1306c write(*,1) m,s,tmn(1) 1307c 1 format('(m,s,tvec)=(',i5,',',i5,',',f15.8,')') 1308c endif 1309 1310 coeff=tmn(1) 1311c ------- get coeff() --- END 1312c Note.- g_munuPSOSO will be the (hyp-para)_uv matrix 1313 jlo=nsize2*(iat1-1)+nsize*(s-1)+1 1314 jhi=jlo+nsize-1 1315 call ga_copy_patch('n',g_munuPSOSO,1,1,jlo,jhi, 1316 & g_tnp1 ,1,1,1 ,nsize) 1317 call ga_add(1.0d0,g_acc1,coeff,g_tnp1,g_acc1) 1318c ----- Calculate rotated perturbed MO: g_acc2 ----- START 1319c \sum_{s=1,3} t_{sm} C_{ri}^{(s) sigma} --> g_acc2 1320 elo(1) = 1 1321 ehi(1) = nbf 1322 elo(2) = 1 1323 ehi(2) = ntot 1324 elo(3) = s 1325 ehi(3) = s 1326 flo(1) = 1 1327 fhi(1) = nbf 1328 flo(2) = 1 1329 fhi(2) = ntot 1330 flo(3) = 3*(iat1-1)+m 1331 fhi(3) = 3*(iat1-1)+m 1332 call nga_add_patch(1.0d0,g_acc2,flo,fhi, 1333 & coeff,g_c1 ,elo,ehi, 1334 & g_acc2,flo,fhi) 1335c ----- Calculate rotated perturbed MO: g_acc2 ----- END 1336 enddo ! end-loop-s 1337c Note: g_acc1 = \widetilde{H}_{mu nu}^{(m)} m=x,y,z 1338c it is the rotated munu matrix using: T H 1339c ====== Store final munu matrices === START 1340 jlo2=nlst*Ndir_munu*(iat1-1)+ 1341 & nlst*(n-1)+1 1342 jhi2=jlo2+nlst-1 1343 call ga_scale(g_acc ,const_hfine(iat1)) ! scaling to conv. to MHz 1344 call ga_scale(g_acc1,const_hfine(iat1)) ! scaling to conv. to MHz 1345 call ga_copy_patch('n',g_acc ,1,1, 1,nlst, 1346 & g_munu_rot,1,1,jlo2,jhi2) 1347 call ga_copy_patch('n',g_acc1 ,1,1, 1,nlst, 1348 & g_munu_rot1,1,1,jlo2,jhi2) 1349c ====== Store final munu matrices === END 1350c ++++++++++++++++++CHECK++++ DIAGONALIZATION ==== START 1351c ==== sum (g_acc .* g_sdens1 + Nuclear CONTRIB) 1352c = TOTAL HYP diagonalized 1353 jlo1=1+nbf 1354 jhi1=nsize 1355 call ga_scale_patch(g_acc,1,1,jlo1,jhi1,2.0d0) 1356 chcdata(n)=ga_ddot(g_acc,g_sdens1) 1357c ++++++++++++++++++CHECK++++ DIAGONALIZATION ==== END 1358c +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1359 enddo ! end-loop-n 1360 if (ga_nodeid().eq.0) then 1361 write(*,23) iat, 1362 & chcdata(1), ! dia-x 1363 & chcdata(2),chcdata(3) ! dia-y,z 1364 23 format(' CHC FCSD(xx,yy,zz)(',i3,')=(', 1365 & f15.8,',',f15.8,',',f15.8,')') 1366 endif 1367c ++++++++++ CHECK diagonalization in para hyperfine ++++++++ START 1368c Note.- Variables defined in zora.fh: 1369c g_CiFull, zora scaling factors 1370c filled out with values in dft_zora_scale 1371c zora switches: do_zora,do_NonRel,not_zora_scale 1372 type_NMR=2 ! =1,2,3=shieldings,hyperfine,gshift 1373 call get_P10_rot( 1374 & g_p10, ! out: Perturbed density matrix (munu nbf x nbf x 3 square matrix) 1375 & type_NMR, ! in: =1,2,3=shieldings,hyperfine,gshift 1376 & g_acc2, ! in: rotated perturbed MO vector 1377 & vectors,g_CiFull, ! in: to build zora scaled MO vector 1378 & iat1, ! in: index for selected atom nr =1,nlist 1379 & nbf,nmo,npol,nocc,nvirt, 1380 & do_zora,do_NonRel,not_zora_scale,rtdb) 1381 call fill_munuPSOSO_1( ! g_munuPSOSO --> g_munuPSOSO2d 1382 & g_munu_rot1 , ! in: array with unique elements 1383 & g_munuPSOSO2d , !out: nbf x nbf x 3 munu matrix for ith atom 1384 & iat1 , ! in: = 1,2,..., nlist 1385 & 2, ! in: type_symm = 1 symm = 2 antisymm 1386 & nbf) 1387c +++++++++ NOW: do ddot product to get diagonalized tensor +++ START 1388 call get_par_hfine_rot( 1389 & g_munuPSOSO2d, ! IN : h01 matrix 1390 & g_p10, ! IN : Perturbed density matrix 1391 & basis,nbf,iat,rtdb) 1392c +++++++++ NOW: do ddot product to get diagonalized tensor +++ END 1393c ++++++++++ CHECK diagonalization in para hyperfine ++++++++ END 1394 if (.not. ga_destroy(g_p10)) call errquit( 1395 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 1396 enddo ! end-loop-iat 1397 if (ga_nodeid().eq.0) 1398 & write(*,*) 'CHCooooooooooooo', 1399 & ' NW-Hyperfine: Summary C+HC data [MHz] ', 1400 & 'oooooooooooooooo END' 1401c --> Main outputs: g_acc2 ,rotated perturbed MO nbf*ntot*ndir*nlist 1402c ndir=1,2,3=x,y,z 1403c nbf, nr of basis functions 1404c ntot=nocc(1)+nocc(2) 1405c nlist, nr of selected atoms 1406c g_munu_rot1,rotated perturbed AO matrix 1407c storing only diag + off-diag elements 1408c Reminder: this comes from an antisymmetrix matrix 1409c in case we want to pull back the 2d munu-matrix 1410c g_munu_rot, rotated AO matrix for dia part 1411c storing only diag + off-diag elements 1412c Reminder: this comes from a symmetric matrix 1413c in case we want to pull back the 2d munu-matrix 1414 call reorder_munu(g_munu_rot ,nat,nlst,nbf,Ndir_munu) ! reoder-munu matrix 1415 call reorder_munu(g_munu_rot1,nat,nlst,nbf,Ndir_munu) ! reoder-munu matrix 1416c ------ destroy unnecessary GAs 1417 if (.not. ga_destroy(g_munuFCSD)) call errquit( 1418 & 'create_munu4nbo-1: ga_destroy failed ',0, GA_ERR) 1419 if (.not. ga_destroy(g_munuPSOSO)) call errquit( 1420 & 'create_munu4nbo-2: ga_destroy failed ',0, GA_ERR) 1421 if (.not. ga_destroy(g_munuPSOSO2d)) call errquit( 1422 & 'create_munu4nbo-7: ga_destroy failed ',0, GA_ERR) 1423 if (.not. ga_destroy(g_tnp)) call errquit( 1424 & 'create_munu4nbo-5: ga_destroy failed ',0, GA_ERR) 1425 if (.not. ga_destroy(g_acc)) call errquit( 1426 & 'create_munu4nbo-6: ga_destroy failed ',0, GA_ERR) 1427 if (.not. ga_destroy(g_tnp1)) call errquit( 1428 & 'create_munu4nbo-5a: ga_destroy failed ',0, GA_ERR) 1429 if (.not. ga_destroy(g_acc1)) call errquit( 1430 & 'create_munu4nbo-6a: ga_destroy failed ',0, GA_ERR) 1431 if (.not. ga_destroy(g_sdens)) call errquit( 1432 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 1433 if (.not. ga_destroy(g_sdens1)) call errquit( 1434 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 1435 if (.not. ga_destroy(g_c1)) call errquit( 1436 & 'create_munu4nbo: ga_destroy failed ',0, GA_ERR) 1437 do i=1,npol 1438 911 if (.not.ga_destroy(vectors(i))) call 1439 & errquit('create_munu4nbo: ga_destroy failed vectors', 1440 & 0,GA_ERR) 1441 enddo 1442c ------ Free allocated memory 1443 call whypfile( 1444 & rtdb, 1445 & g_munu_rot, 1446 & g_munu_rot1, 1447 & g_acc2, 1448 & nlst, 1449 & Ndir_munu, 1450 & Natoms_munu, 1451 & atmnr_munu) 1452 if (.not. ga_destroy(g_munu_rot)) call errquit( 1453 & 'whypfile: ga_destroy failed ',0, GA_ERR) 1454 if (.not. ga_destroy(g_munu_rot1)) call errquit( 1455 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 1456 if (.not. ga_destroy(g_acc2)) call errquit( 1457 & 'wshldfile: ga_destroy failed ',0, GA_ERR) 1458 return 1459 end 1460 1461 subroutine get_P10_rot( 1462 & g_p1, ! out: Perturbed (spin)density matrix 1463 & type_NMR, ! in: =1,2,3=shieldings,hyperfine,gshift 1464 & g_c1, ! in: rotated perturbed MO 1465 & vectors, ! in: MO coeffs 1466 & g_CiFull, ! in: MO zora weighting coeffs 1467 & iat1, ! in: index of selected atoms =1,nlist 1468 & nbf, ! in: nr. basis functions 1469 & nmo, ! in: nr. MOs (occ+virt) 1470 & npol, ! in: nr. of polarizations 1471 & nocc, ! in: nr. occ MOs 1472 & nvirt, ! in: nr. virtual MOs 1473 & do_zora, ! in: .true. if doing zora calc 1474 & do_NonRel, ! in: .true. if doing nonrel within zora scheme 1475 & not_zora_scale, ! in: .true. not scaling perturbed density matrix 1476 & rtdb) 1477c Note: If shieldings calc --> g_p1= Perturbed density matrix 1478c = P^(1,0)_A + P^(1,0)_B 1479c ==> This calc. could be npol=1 ( restricted shell calc.) 1480c npol=2 (unrestricted shell calc.) 1481c If gshift calc --> g_p1= Perturbed (spin) density matrix 1482c = P^(1,0)_A - P^(1,0)_B 1483c Purpose: Compute rotated P10 (perturbed density matrix) 1484c Using a rotated perturbed MO (g_c1) 1485c This routine is modified from get_P10 1486c and is intended to be used in hyperfine coupling constants 1487c NLMO analysis. 1488c Author : Fredy W. Aquino 1489c Date : 07-07-11 1490 implicit none 1491#include "errquit.fh" 1492#include "global.fh" 1493#include "mafdecls.fh" 1494#include "msgids.fh" 1495#include "geom.fh" 1496#include "rtdb.fh" 1497#include "bas.fh" 1498#include "stdio.fh" 1499 integer npol ! nr. of polarizations 1500 integer g_p1 ! OUT: Perturbed density matrix 1501 integer g_u,g_d1 ! scratch ga array 1502 integer vectors(npol),vectors_scl(npol), 1503 & g_CiFull(npol) 1504 integer rtdb 1505 integer nbf,nmo,ispin,type_NMR, 1506 & nocc(npol),nvirt(npol) 1507 integer shift,xyz,iat1 1508 integer alo(3), ahi(3), 1509 & blo(3), bhi(3), 1510 & clo(3), chi(3), 1511 & dlo(3), dhi(3), 1512 & elo(3), ehi(3) 1513 double precision coeff(2) 1514 logical status,do_zora,do_NonRel,not_zora_scale 1515c ------ for Hyperfine NMLO analysis --------- START 1516 integer ntot,g_c1, ! g_c1 , collects perturbed MO coeffs C1 1517 & plo(3),phi(3),qlo(3),qhi(3), 1518 & nlist 1519c ------ for Hyperfine NMLO analysis --------- END 1520 if (type_NMR.eq.1) then ! Shieldings 1521 coeff(1)= 1.0d0 1522 coeff(2)= 1.0d0 1523 else if (type_NMR.eq.2 .or. ! Hyperfine 1524 & type_NMR.eq.3) then ! g-shifts 1525 coeff(1)= 1.0d0 1526 coeff(2)= -1.0d0 1527 else 1528 write(*,*) 'Error in get_P10_1:', 1529 & ' Calc. should be giao, gshift or hyperfine.' 1530 stop 1531 endif 1532 alo(1) = nbf 1533 alo(2) = -1 1534 alo(3) = -1 1535 ahi(1) = nbf 1536 ahi(2) = nbf 1537 ahi(3) = 3 1538 if (.not.nga_create(MT_DBL,3,ahi,'get_P10_rot: g_d1 matrix', 1539 & alo,g_d1)) call 1540 & errquit('g_d1: nga_create failed g_d1',0,GA_ERR) 1541 if (.not.nga_create(MT_DBL,3,ahi,'get_P10_rot: g_p1 matrix', 1542 & alo,g_p1)) call 1543 & errquit('g_p1: nga_create failed g_p1',0,GA_ERR) 1544 call ga_zero(g_p1) 1545 do ispin=1,npol 1546 alo(1) = nbf 1547 alo(2) = -1 1548 alo(3) = -1 1549 ahi(1) = nbf 1550 ahi(2) = nocc(ispin) 1551 ahi(3) = 3 1552 if (.not.nga_create(MT_DBL,3,ahi,'U matrix',alo,g_u)) call 1553 & errquit('g_u: nga_create failed g_u',0,GA_ERR) 1554 call ga_zero(g_u) 1555c From U matrices, generate the perturbed density matrices D1x,y,z 1556c C1 = C0 * U10 1557c D1 = 2[(C1*C0+) - (C0*C1+)] 1558 alo(1) = 1 1559 alo(2) = 1 1560 blo(1) = 1 1561 blo(2) = 1 1562 clo(1) = 1 1563 chi(1) = nbf 1564 clo(2) = 1 1565 chi(2) = nbf 1566 dlo(1) = 1 1567 dlo(2) = 1 1568 dhi(1) = nbf 1569 dhi(2) = nocc(ispin) 1570c --------- zora scaling of MO vectors(1) ----- START 1571c Note.- g_CiFull is defined in dft_zora_scale() (source dft_zora_utils.F) 1572 if(.not.ga_duplicate(vectors(ispin), 1573 & vectors_scl(ispin),'vscl 1')) 1574 & call errquit('g_d1: ga_duplicate failed',1,GA_ERR) 1575 call ga_copy(vectors(ispin),vectors_scl(ispin)) 1576 if (do_zora .and. .not.(do_NonRel) .and. 1577 & .not.(not_zora_scale)) then 1578 call ga_scale_cols(vectors_scl(ispin),g_CiFull(ispin)) 1579 endif 1580c --------- zora scaling of MO vectors(1) ----- END 1581 do xyz = 1,3 ! = x,y,z 1582 alo(3) = xyz 1583 ahi(3) = xyz 1584 blo(3) = xyz 1585 bhi(3) = xyz 1586 clo(3) = xyz 1587 chi(3) = xyz 1588 dlo(3) = xyz 1589 dhi(3) = xyz 1590 bhi(1) = nbf 1591 bhi(2) = nmo 1592 ahi(1) = nmo 1593 ahi(2) = nocc(ispin) 1594c Make C1 1595c +++++++++++ here g_c1 --> g_u +++++++++++ START 1596 shift=nocc(1)*(ispin-1) 1597 qlo(1) = 1 1598 qhi(1) = nbf 1599 qlo(2) = shift+1 1600 qhi(2) = shift+nocc(ispin) 1601 qlo(3) = 3*(iat1-1)+xyz 1602 qhi(3) = 3*(iat1-1)+xyz 1603 plo(1) = 1 1604 phi(1) = nbf 1605 plo(2) = 1 1606 phi(2) = nocc(ispin) 1607 plo(3) = xyz 1608 phi(3) = xyz 1609 call nga_copy_patch('n',g_c1,qlo,qhi, 1610 & g_u ,plo,phi) 1611c +++++++++++ here g_c1 --> g_u +++++++++++ END 1612 bhi(1) = nbf 1613 bhi(2) = nocc(ispin) 1614 ahi(1) = nocc(ispin) 1615 ahi(2) = nbf 1616c Make D1 1617 call nga_matmul_patch('n','t',1.0d0,0.0d0, 1618 & vectors_scl(ispin),blo,bhi, 1619 & g_u,alo,ahi, 1620 & g_d1,clo,chi) 1621 call nga_matmul_patch('n','t',-1.0d0,1.0d0, 1622 & g_u,blo,bhi, 1623 & vectors_scl(ispin),alo,ahi, 1624 & g_d1,clo,chi) 1625 enddo ! end-loop-xyz 1626 elo(1) = 1 1627 ehi(1) = nbf 1628 elo(2) = 1 1629 ehi(2) = nbf 1630 elo(3) = 1 1631 ehi(3) = 3 1632 call nga_add_patch(1.0d0 ,g_p1,elo,ehi, 1633 & coeff(ispin),g_d1,elo,ehi, 1634 & g_p1,elo,ehi) 1635 if (.not.ga_destroy(g_u)) call 1636 & errquit('get_d1: ga_destroy failed g_u',0,GA_ERR) 1637 if (.not.ga_destroy(vectors_scl(ispin))) call 1638 & errquit('get_d1: ga_destroy failed vscl',0,GA_ERR) 1639 enddo ! end-loop-ispin 1640 if (npol.eq.1 .and. type_NMR.eq.1) then ! this happens ONLY for NMR-restricted calc. 1641 elo(1) = 1 1642 ehi(1) = nbf 1643 elo(2) = 1 1644 ehi(2) = nbf 1645 elo(3) = 1 1646 ehi(3) = 3 1647 call nga_scale_patch(g_p1,elo,ehi,2.0d0) 1648 endif 1649 if (.not.ga_destroy(g_d1)) call 1650 & errquit('get_d1: ga_destroy failed g_d1',0,GA_ERR) 1651 return 1652 end 1653 1654 subroutine get_par_hfine_rot( ! for checking rotated Pert. MO and rotated munu AO 1655 & ga_h01_hfine,! IN : 1656 & g_d1, ! IN : Perturbed density matrix 1657 & basis, ! IN : basis handle 1658 & nbf, ! IN : nr. basis functions 1659 & iat, ! IN : atom nr 1660 & rtdb) ! IN : rtdb handle 1661c Purpose : Assemble NMR Hyperfine: paramagnetic tensor 1662c Author : Fredy Aquino 1663c Date : 07-08-11 1664c Note : This routine is taken from get_par_hfine() 1665c The only difference is that here natoms=1 1666 implicit none 1667#include "errquit.fh" 1668#include "global.fh" 1669#include "mafdecls.fh" 1670#include "msgids.fh" 1671#include "rtdb.fh" 1672#include "apiP.fh" 1673#include "prop.fh" 1674#include "bgj.fh" 1675#include "bas.fh" 1676#include "stdio.fh" 1677 1678 integer g_d1 ! input: Perturbed density matrices 1679 integer ga_h01_hfine 1680 integer rtdb 1681 integer basis 1682 integer alo(3), ahi(3) 1683 integer ld(2),cbuf 1684 integer nbf,iat 1685 integer iy 1686 double precision val,chcdata(3) 1687 logical oskel 1688 1689 alo(1) = 1 1690 ahi(1) = nbf 1691 alo(2) = 1 1692 ahi(2) = nbf 1693 do iy = 1, 3 1694 alo(3) = iy 1695 ahi(3) = iy 1696 val = nga_ddot_patch(g_d1 ,'n',alo,ahi, 1697 & ga_h01_hfine,'n',alo,ahi) 1698 chcdata(iy)=val 1699 enddo ! end-loop-iy 1700 if (ga_nodeid().eq.0) then 1701 write(*,1) iat, 1702 & chcdata(1), ! par-x 1703 & chcdata(2),chcdata(3) ! par-y,z 1704 1 format(' CHC PSOSO(xx,yy,zz)(',i3,')=(', 1705 & f15.8,',',f15.8,',',f15.8,')') 1706 endif 1707 return 1708 end 1709 subroutine get2dmat(a,n,b) 1710c a : size(a)=n*(n+1)/2) symmetric matrix 1711 implicit none 1712#include "nwc_const.fh" 1713#include "errquit.fh" 1714#include "global.fh" 1715#include "bas.fh" 1716#include "mafdecls.fh" 1717#include "geom.fh" 1718#include "stdio.fh" 1719#include "msgids.fh" 1720 integer a,b,n,n1,i,j,count, 1721 & l_mat,k_mat 1722 1723 if (.not. ga_create(mt_dbl,n,n, 1724 & 'get2dmat: b',0,0,b)) 1725 $ call errquit('get2dmat: b', 0,GA_ERR) 1726 call ga_zero(b) 1727 1728 n1=n*(n+1)/2 1729 if (.not.ma_alloc_get(mt_dbl,n1, 1730 & 'fcsd',l_mat,k_mat)) 1731 & call errquit('get2dmat: ma failed',0,MA_ERR) 1732 call ga_get(a,1,1,1,n1, 1733 & dbl_mb(k_mat),1) 1734 count=0 1735 do i=1,n 1736 call ga_put(b,i,i,i,i,dbl_mb(k_mat+count),1) 1737 count=count+1 1738 enddo 1739 do i=2,n 1740 do j=1,i-1 1741 call ga_put(b,i,i,j,j,dbl_mb(k_mat+count),1) 1742 call ga_put(b,j,j,i,i,dbl_mb(k_mat+count),1) 1743 count=count+1 1744 enddo ! end-loop-j 1745 enddo ! end-loop-i 1746 if (.not.ma_free_heap(l_mat)) call 1747 & errquit('get2dmat:: ma_free_heap l_mat',0, MA_ERR) 1748 return 1749 end 1750c -------------- for NMLO analysis ----------------- END 1751c -------------------------------------------------------- 1752