1 2*********************************************************************** 3* cpsdv5-mpi (MPI code) * 4* * 5* This is a developing cpsdv5 parallel code for the SP2 * 6* + mpi message passing library used * 7* + ngp is used instead of nfft in this proceudure * 8* + error checking is based on aimd.h parameters * 9* then control file * 10* + my own slap-decomposed parallel 3d-FFT(real->complex) used * 11* * 12* * 13*********************************************************************** 14 15 logical function cpsdv5(rtdb) 16* 17* $Id$ 18* 19 implicit none 20 integer rtdb 21 22#include "global.fh" 23#include "bafdecls.fh" 24#include "stdio.fh" 25#include "util.fh" 26#include "btdb.fh" 27#include "errquit.fh" 28cccc#include "frac_occ.fh" 29 30* **** parallel variables **** 31 logical omaster 32 integer taskid,np,np_i,np_j 33 integer MASTER 34 parameter(MASTER=0) 35 36* **** timing variables **** 37 real*8 cpu1,cpu2,cpu3,cpu4 38 real*8 t1,t2,t3,t4,av 39 40* **** lattice variables **** 41 integer ngrid(3),nwave,nfft3d,n2ft3d,ngrid_small(3) 42 integer npack1 43 real*8 a,b,c,alpha,beta,gamma 44 45 46* **** electronic variables **** 47 logical first_iteration 48 real*8 icharge 49 integer ispin 50 integer ne(2),n1(2),n2(2),nemax,neq(2),nemaxq 51 real*8 en(2),en1(2),en2(2) 52 real*8 dipole(3) 53 54* complex*16 psi1(nfft3d,nemax) 55* complex*16 psi2(nfft3d,nemax) 56* real*8 dn(n2ft3d,2) 57* complex*16 Hpsi(nfft3d,nemax) 58* complex*16 psir(nfft3d,nemax) 59 integer psi1(2),psi2(2) 60 integer occ1(2),occ2(2) 61 integer dn(2) 62 integer Hpsi(2),psir(2) 63 64 65* ***** energy variables **** 66 real*8 E(60),Egas 67 68* real*8 eig(2*nemax) 69* real*8 hml(2*nemax*nemax) 70* real*8 lmd(2*nemax*nemax) 71 integer eig(2),hml(2),lmd(2),lmd1(2) 72 73* **** psi smearing block **** 74 logical fractional 75 integer smearoccupation,smeartype 76 real*8 smearfermi(2),smearcorrection,smearkT 77 78 79 80 81* **** error variables **** 82 integer ierr 83 84* **** local variables **** 85 logical oprint,lprint,hprint,mparallelized 86 integer ms,mapping,mapping1d 87 real*8 deltae,deltac,deltar 88 real*8 gx,gy,gz,cx,cy,cz,sum1,sum2 89 real*8 EV,pi 90 integer i,j,k,ia,n,nn 91 integer ii,jj,indx,indx1 92 integer icount,it_in,it_out 93 real*8 w,sumall,virial,dv 94 integer nfft3 95 parameter (nfft3=32) 96 character*255 full_filename 97 integer tmp1(2) 98 99 logical value,psi_nogrid,field_exist 100 integer hversion,hnfft(3),hispin,hne(2) 101 real*8 hunita(3,3) 102 integer ind 103 character*50 filename 104 character*50 control_input_psi 105 external control_input_psi 106 logical wvfnc_expander 107 external wvfnc_expander 108 109 110 111 112* **** external functions **** 113 real*8 psp_zv,psp_rc,ewald_rcut,ion_amass,nwpw_cosmo_qc0 114 real*8 ewald_mandelung,lattice_omega_small 115 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 116 real*8 lattice_unitg,lattice_unitg_small,lattice_unita_small 117 integer ewald_ncut,ewald_nshl3d 118 integer psp_lmmax,psp_lmax,psp_locp 119 character*4 ion_aname,ion_atom 120 external psp_zv,psp_rc,ewald_rcut,ion_amass,nwpw_cosmo_qc0 121 external ewald_mandelung,lattice_omega_small 122 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 123 external lattice_unitg,lattice_unitg_small,lattice_unita_small 124 external ewald_ncut,ewald_nshl3d 125 external psp_lmmax,psp_lmax,psp_locp 126 external ion_aname,ion_atom 127 128 character*50 control_cell_name 129 external control_cell_name 130 131 real*8 control_tole,control_tolc,control_tolr,ion_rion 132 external control_tole,control_tolc,control_tolr,ion_rion 133 real*8 control_time_step,control_fake_mass 134 external control_time_step,control_fake_mass 135 logical control_read,control_move,ion_init,ion_q_FixIon 136 external control_read,control_move,ion_init,ion_q_FixIon 137 logical ion_q_xyzFixIon 138 external ion_q_xyzFixIon 139 character*14 ion_q_xyzFixIon_label 140 external ion_q_xyzFixIon_label 141 142 integer pack_nwave_all 143 integer control_it_in,control_it_out,control_gga,control_version 144 integer control_ngrid,pack_nwave 145 integer ion_nion,ion_natm,ion_katm,ion_nkatm,ion_nkatm0 146 integer ion_nion_qm,ion_nion_mm,ion_nion_q 147 external pack_nwave_all 148 external control_it_in,control_it_out,control_gga,control_version 149 external control_ngrid,pack_nwave 150 external ion_nion,ion_natm,ion_katm,ion_nkatm,ion_nkatm0 151 external ion_nion_qm,ion_nion_mm,ion_nion_q 152 153 character*12 control_boundry 154 external control_boundry 155 156 logical pspw_HFX_localize2,pspw_charge_found,pspw_Efield_found 157 external pspw_HFX_localize2,pspw_charge_found,pspw_Efield_found 158 logical pspw_SIC,pspw_SIC_relaxed,pspw_qmmm_found 159 logical pspw_HFX,pspw_HFX_relaxed,psp_pawexist 160 logical psp_semicore,control_Mulliken 161 real*8 psp_rcore,psp_ncore,psp_rlocal 162 external pspw_SIC,pspw_SIC_relaxed,pspw_qmmm_found 163 external pspw_HFX,pspw_HFX_relaxed,psp_pawexist 164 external psp_semicore,control_Mulliken 165 external psp_rcore,psp_ncore,psp_rlocal 166 logical control_check_charge_multiplicity 167 external control_check_charge_multiplicity 168 real*8 nwpw_timing 169 external nwpw_timing 170 integer control_np_orbital,control_mapping,control_mapping1d 171 external control_np_orbital,control_mapping,control_mapping1d 172 173 logical control_translation,control_rotation,control_balance 174 external control_translation,control_rotation,control_balance 175 176 logical Dneall_m_allocate,Dneall_m_free,control_parallel_io 177 external Dneall_m_allocate,Dneall_m_free,control_parallel_io 178 179 real*8 Dneall_m_value 180 external Dneall_m_value 181 character*9 ion_amm 182 external ion_amm 183 logical psp_U_psputerm,meta_found,ion_disp_on 184 external psp_U_psputerm,meta_found,ion_disp_on 185 integer ion_nconstraints,ion_ndof 186 external ion_nconstraints,ion_ndof 187 logical control_print,ion_makehmass2,control_has_ngrid_small 188 external control_print,ion_makehmass2,control_has_ngrid_small 189 integer control_ngrid_small 190 external control_ngrid_small 191c#include "Parallel.fh" 192 integer Parallel_threadid,Parallel_nthreads,Parallel_maxthreads 193 external Parallel_threadid,Parallel_nthreads,Parallel_maxthreads 194 integer tid 195 196 logical control_mparallelized 197 external control_mparallelized 198 logical nwpw_cosmo_on,nwpw_born_on 199 external nwpw_cosmo_on,nwpw_born_on 200 real*8 nwpw_cosmo_screen,nwpw_born_screen 201 external nwpw_cosmo_screen,nwpw_born_screen 202 logical pspw_V_APC_on 203 external pspw_V_APC_on 204 real*8 control_gas_energy 205 external control_gas_energy 206 207 logical control_fractional 208 external control_fractional 209 integer control_fractional_smeartype 210 external control_fractional_smeartype 211 real*8 control_fractional_temperature 212 external control_fractional_temperature 213 real*8 control_fractional_kT 214 external control_fractional_kT 215 real*8 control_fractional_alpha 216 external control_fractional_alpha 217 218 integer pspw_Efield_type 219 external pspw_Efield_type 220 221 222 223* |************| 224*****************************| PROLOGUE |**************************** 225* |************| 226 227 value = .true. 228 pi = 4.0d0*datan(1.0d0) 229 230 call nwpw_timing_init() 231 call dcopy(60,0.0d0,0,E,1) 232 233 234* **** get parallel variables **** 235 call Parallel_Init() 236 call Parallel_np(np) 237 call Parallel_taskid(taskid) 238 239 value = MA_set_hard_fail(.true.) 240 241 tid = Parallel_threadid() 242 omaster = (taskid.eq.MASTER).and.(Parallel_threadid().eq.MASTER) 243 244 if (omaster) call current_second(cpu1) 245 if (omaster) call current_second(cpu1) 246 247* ***** print out header **** 248 if (omaster) then 249 write(luout,1000) 250 write(luout,1010) 251 write(luout,1020) 252 write(luout,1010) 253 write(luout,1030) 254 write(luout,1010) 255 write(luout,1035) 256 write(luout,1010) 257 write(luout,1040) 258 write(luout,1010) 259 write(luout,1041) 260 write(luout,1042) 261 write(luout,1043) 262 write(luout,1010) 263 write(luout,1000) 264 call nwpw_message(1) 265 write(luout,1110) 266 end if 267 268 value = control_read(1,rtdb) 269 270 271 mparallelized = control_mparallelized() 272 call Parallel2d_Init(control_np_orbital()) 273 call Parallel2d_np_i(np_i) 274 call Parallel2d_np_j(np_j) 275 276 ngrid(1) = control_ngrid(1) 277 ngrid(2) = control_ngrid(2) 278 ngrid(3) = control_ngrid(3) 279 nwave = 0 280 mapping = control_mapping() 281 282 oprint = omaster.and.control_print(print_medium) 283 lprint = omaster.and.control_print(print_low) 284 hprint = omaster.and.control_print(print_high) 285 286 287* **** initialize psi_data **** 288 call psi_data_init(100) 289 290* **** initialize D3dB data structure **** 291 call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 292 call D3dB_nfft3d(1,nfft3d) 293 n2ft3d = 2*nfft3d 294 295* ***** Initialize double D3dB data structure **** 296 if (control_version().eq.4) 297 > call D3dB_Init(2,2*ngrid(1),2*ngrid(2),2*ngrid(3),mapping) 298 299 if (control_has_ngrid_small()) then 300 ngrid_small(1) = control_ngrid_small(1) 301 ngrid_small(2) = control_ngrid_small(2) 302 ngrid_small(3) = control_ngrid_small(3) 303 call D3dB_Init(3,ngrid_small(1),ngrid_small(2),ngrid_small(3), 304 > mapping) 305 end if 306 307 308* **** initialize lattice and packing data structure **** 309 call lattice_init() 310 call G_init() 311 call mask_init() 312 call Pack_Init() 313 call Pack_npack(1,npack1) 314 call D3dB_pfft_init() 315 !call ga_sync() 316 317* **** read ions **** 318 value = ion_init(rtdb) 319 call center_geom(cx,cy,cz) 320 call center_mass(gx,gy,gz) 321 first_iteration = .true. 322 323* **** allocate psp data structure and read in psedupotentials into it **** 324 call psp_init() 325 call psp_readall() 326 if (psp_semicore(0)) call semicore_check() 327 328* **** initialize G,mask,ke,and coulomb data structures **** 329 call ke_init() 330 if (control_version().eq.3) call coulomb_init() 331 if (control_version().eq.4) call coulomb2_init() 332 call strfac_init() 333 call phafac() 334 if (control_version().eq.3) call ewald_init() 335 336* **** generate initial wavefunction if it does not exist **** 337 if (.not.control_check_charge_multiplicity()) then 338 call psi_new() 339 end if 340 341* ***** allocate psi2,and psi1 wavefunctions **** 342 call psi_get_ne_occupation(ispin,ne,smearoccupation) 343 if (smearoccupation.gt.0) then 344 fractional = .true. 345 else 346 fractional = .false. 347 end if 348 smearcorrection = 0.0d0 349 smeartype = control_fractional_smeartype() 350 smearkT = control_fractional_kT() 351 352 mapping1d = control_mapping1d() 353 call Dne_init(ispin,ne,mapping1d) 354 call Dneall_neq(neq) 355 nemaxq = neq(1)+neq(2) 356 357 value = BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 358 > 'psi2',psi2(2),psi2(1)) 359 value = value.and. 360 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 361 > 'psi1',psi1(2),psi1(1)) 362 if (fractional) then 363 value = value.and. 364 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ1',occ1(2),occ1(1)) 365 value = value.and. 366 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ2',occ2(2),occ2(1)) 367 end if 368 if (.not. value) call errquit('out of heap memory',0, MA_ERR) 369 370 371* ***** read initial wavefunctions into psi2 **** 372 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 373 > mt_log,1,psi_nogrid)) 374 > psi_nogrid = .true. 375 376 if (psi_nogrid) then 377 378 call psi_get_header(hversion,hnfft,hunita,hispin,hne) 379 380 if ( (hnfft(1).ne.control_ngrid(1)) .or. 381 > (hnfft(2).ne.control_ngrid(2)) .or. 382 > (hnfft(3).ne.control_ngrid(3)) ) then 383 384 hnfft(1) = control_ngrid(1) 385 hnfft(2) = control_ngrid(2) 386 hnfft(3) = control_ngrid(3) 387 call Parallel_taskid(taskid) 388 389 call ga_sync() 390 value = btdb_parallel(.false.) 391 call ga_sync() 392 if (oprint) then 393 394 filename = control_input_psi() 395 396 ind = index(filename,' ') - 1 397 if (.not. btdb_cput(rtdb,'xpndr:old_wavefunction_filename', 398 > 1,filename(1:ind))) 399 > call errquit( 400 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 401 402 if (.not. btdb_cput(rtdb,'xpndr:new_wavefunction_filename', 403 > 1,filename(1:ind))) 404 > call errquit( 405 > 'wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 406 407 if (.not. btdb_put(rtdb,'xpndr:ngrid',mt_int,3,hnfft)) 408 > call errquit( 409 > 'wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 410 411 write(*,*) 412 write(*,*) "Grid is being converted:" 413 write(*,*) "------------------------" 414 write(*,*) 415 write(*,*) "To turn off automatic grid conversion:" 416 write(*,*) 417 write(*,*) "set nwpw:psi_nogrid .false." 418 write(*,*) 419 value = wvfnc_expander(rtdb) 420 421 end if 422 call ga_sync() 423 value = btdb_parallel(.true.) 424 425 end if 426 427 end if 428 429* ***** read psi2 wavefunctions **** 430 call psi_read(ispin,ne,dcpl_mb(psi2(1)), 431 > smearoccupation,dbl_mb(occ2(1))) 432 if (smearoccupation.gt.0) 433 > call Parallel_shared_vector_copy(.true.,ne(1)+ne(2), 434 > dbl_mb(occ2(1)), 435 > dbl_mb(occ1(1))) 436 437 438 n1(1) = 1 439 n2(1) = ne(1) 440 n1(2) = ne(1)+1 441 n2(2) = ne(1)+ne(2) 442 nemax = ne(1)+ne(2) 443 444* **** Ortho Check **** 445 if (psp_pawexist()) then 446 call psp_overlap_S(ispin,neq,dcpl_mb(psi2(1)),dcpl_mb(psi1(1))) 447 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 448 > dcpl_mb(psi1(1)), 449 > dcpl_mb(psi2(1)), 450 > sum2) 451 else 452 call Grsm_gg_trace(npack1,(neq(1)+neq(2)), 453 > dcpl_mb(psi2(1)), 454 > dcpl_mb(psi2(1)), 455 > sum2) 456 end if 457 call D1dB_SumAll(sum2) 458 459 sum1 = dble(ne(1) + ne(2)) 460 if (dabs(sum2-sum1).gt.1.0d-10) then 461 if (oprint) 462 > write(*,*) "Warning: Gram-Schmidt being performed on psi2" 463 if (psp_pawexist()) then 464 call Dneall_f_Sortho(0,dcpl_mb(psi2(1)), 465 > dcpl_mb(psi1(1)),npack1) 466 else 467 call Dneall_f_GramSchmidt(0,dcpl_mb(psi2(1)),npack1) 468 end if 469c call Dneall_f_ortho(0,dcpl_mb(psi2(1)),npack1) 470c call Grsm_g_MakeOrtho(npack1,ne(1),dcpl_mb(psi2(1))) 471c if (ispin.gt.1) then 472c call Grsm_g_MakeOrtho(npack1,ne(2), 473c > dcpl_mb(psi2(1) + ne(1)*npack1)) 474c end if 475 end if 476 477* **** allocate other variables ***** 478 value = BA_alloc_get(mt_dbl,(2*nemax),'eig',eig(2),eig(1)) 479 value = value.and.Dneall_m_allocate(0,hml) 480 value = value.and.Dneall_m_allocate(0,lmd) 481 value = value.and.Dneall_m_allocate(0,lmd1) 482 483 value = value.and. 484 > BA_alloc_get(mt_dbl,(4*nfft3d), 485 > 'dn',dn(2),dn(1)) 486 value = value.and. 487 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 488 > 'Hpsi',Hpsi(2),Hpsi(1)) 489 value = value.and. 490 > BA_alloc_get(mt_dcpl,nfft3d*(neq(1)+neq(2)), 491 > 'psir',psir(2),psir(1)) 492 493 if (.not. value) 494 > call errquit('cpsdv5:out of heap memory',0,MA_ERR) 495 496 497 498* **** initialize two-electron Gaussian integrals **** 499* **** initialize paw ncmp*Vloc **** 500 if (psp_pawexist()) then 501 call nwpw_gintegrals_init() 502 call nwpw_gintegrals_set(control_move()) 503 call psp_dE_ncmp_vloc_Qlm(ispin,.false.,hunita) 504 end if 505 506 507* **** initialize SIC and HFX **** 508 call pspw_init_SIC(rtdb,ne) 509 call pspw_init_HFX(rtdb,ispin,ne) 510 511* **** initialize DFT+U **** 512 call psp_U_init() 513 514* **** initialize META GGA **** 515 call nwpw_meta_gga_init(control_gga()) 516 517* **** initialize vdw **** 518 call vdw_DF_init() 519 520* **** initialize rho_symmetry **** 521 call rho_symmetrizer_init() 522 523* **** initialize metadynamics **** 524 call meta_initialize(rtdb) 525 526* **** initialize QM/MM **** 527 call pspw_init_APC(rtdb) 528 call pspw_qmmm_init(rtdb) 529 call pspw_charge_init(rtdb) 530 call pspw_Efield_init(rtdb,ispin,ne) 531 !call pspw_cdft_init(rtdb) 532 field_exist = pspw_charge_found().or.pspw_Efield_found() 533 534* **** initialize frac_occ data structure **** 535c call frac_occ_init(rtdb,ispin,ne) 536 537* **** initialize FixIon constraint **** 538 call ion_init_FixIon(rtdb) 539 540 541 542* |**************************| 543****************** summary of input data ********************** 544* |**************************| 545 546 if (oprint) then 547 write(luout,1111) np 548 write(luout,1117) np_i,np_j 549 if (mapping.eq.1) write(luout,1112) 550 if (mapping.eq.2) write(luout,1113) 551 if (mapping.eq.3) write(luout,1118) 552 if (control_balance()) then 553 write(luout,1114) 554 else 555 write(luout,1116) 556 end if 557 if (control_parallel_io()) then 558 write(luout,1119) 559 else 560 write(luout,1122) 561 end if 562 write(luout,1123) Parallel_maxthreads() 563 564 565 write(luout,1115) 566 IF(control_move()) THEN 567 write(luout,1120) 'yes' 568 ELSE 569 write(luout,1120) 'no' 570 ENDIF 571 write(luout,1121) control_boundry(),control_version() 572 if (ispin.eq.1) write(luout,1130) 'restricted' 573 if (ispin.eq.2) write(luout,1130) 'unrestricted' 574 575 call v_bwexc_print(luout,control_gga()) 576 577 if (fractional) write(luout,1132) 578 call pspw_print_SIC(luout) 579 call pspw_print_HFX(luout) 580 call nwpw_meta_gga_print(luout) 581 if (ion_makehmass2()) write(luout,1135) 582 write(luout,1140) 583 do ia = 1,ion_nkatm() 584 call psp_print(ia) 585c write(luout,1150) ia,ion_atom(ia), 586c > psp_zv(ia),psp_lmax(ia) 587c write(luout,1152) psp_lmax(ia) 588c write(luout,1153) psp_locp(ia) 589c write(luout,1154) psp_lmmax(ia) 590c if (control_version().eq.4) write(luout,1156) psp_rlocal(ia) 591c if (psp_semicore(ia)) 592c > write(luout,1155) psp_rcore(ia),psp_ncore(ia) 593c write(luout,1151) (psp_rc(i,ia),i=0,psp_lmax(ia)) 594 end do 595 596 icharge = -(ne(1)+ne(ispin)) 597 en(1) = ne(1) 598 en(ispin) = ne(ispin) 599 if (fractional) then 600 icharge = 0.0d0 601 do ms=1,ispin 602 en(ms) =0.0 603 do i=n1(ms),n2(ms) 604 icharge = icharge - (3-ispin)*dbl_mb(occ1(1)+i-1) 605 en(ms) = en(ms) + dbl_mb(occ1(1)+i-1) 606 end do 607 end do 608 end if 609 610 do ia=1,ion_nkatm0() 611 icharge = icharge + ion_natm(ia)*psp_zv(ia) 612 end do 613 write(luout,1159) icharge 614 615 write(luout,1160) 616 write(luout,1170) (ion_atom(K),ion_natm(K),K=1,ion_nkatm()) 617 write(luout,1180) 618 do I=1,ion_nion()+ion_nion_q() 619 if (I.gt.ion_nion()) then 620 write(luout,1195) I,ion_aname(I),(ion_rion(K,I),K=1,3), 621 > nwpw_cosmo_qc0(I-ion_nion()),ion_amm(i) 622 else if (ion_q_FixIon(I)) then 623 write(luout,1191) I,ion_aname(I),(ion_rion(K,I),K=1,3), 624 > ion_amass(I)/1822.89d0,ion_amm(i) 625 else if (ion_q_xyzFixIon(I)) then 626 write(luout,1194) I,ion_aname(I),(ion_rion(K,I),K=1,3), 627 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 628 else 629 write(luout,1190) I,ion_aname(I),(ion_rion(K,I),K=1,3), 630 > ion_amass(I)/1822.89d0,ion_amm(i) 631 end if 632 end do 633 write(luout,1200) cx,cy,cz 634 write(luout,1210) gx,gy,gz 635 write(luout,1211) ion_nconstraints(),ion_ndof() 636 637 call pspw_charge_Print(luout) 638 call pspw_Efield_Print(luout) 639 640 if (fractional) then 641 write(luout,1219) en(1),en(ispin),' ( fractional)' 642 write(luout,1221) ne(1),neq(1), 643 > ne(ispin),neq(ispin),' (Fourier space)' 644 else 645 write(luout,1220) ne(1),neq(1), 646 > ne(ispin),neq(ispin),' (Fourier space)' 647 write(luout,1221) ne(1),neq(1), 648 > ne(ispin),neq(ispin),' (Fourier space)' 649 end if 650 write(luout,1230) 651 write(luout,1241) lattice_unita(1,1), 652 > lattice_unita(2,1), 653 > lattice_unita(3,1) 654 write(luout,1242) lattice_unita(1,2), 655 > lattice_unita(2,2), 656 > lattice_unita(3,2) 657 write(luout,1243) lattice_unita(1,3), 658 > lattice_unita(2,3), 659 > lattice_unita(3,3) 660 write(luout,1244) lattice_unitg(1,1), 661 > lattice_unitg(2,1), 662 > lattice_unitg(3,1) 663 write(luout,1245) lattice_unitg(1,2), 664 > lattice_unitg(2,2), 665 > lattice_unitg(3,2) 666 write(luout,1246) lattice_unitg(1,3), 667 > lattice_unitg(2,3), 668 > lattice_unitg(3,3) 669 write(luout,1231) lattice_omega() 670 write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 671 > pack_nwave_all(0),pack_nwave(0) 672 write(luout,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 673 > pack_nwave_all(1),pack_nwave(1) 674 if (control_version().eq.3) then 675 write(luout,1260) ewald_rcut(),ewald_ncut() 676 write(luout,1261) ewald_mandelung() 677 end if 678 679 if (control_has_ngrid_small()) then 680 write(luout,1229) 681 write(luout,1233) control_cell_name() 682 write(luout,1241) lattice_unita_small(1,1), 683 > lattice_unita_small(2,1), 684 > lattice_unita_small(3,1) 685 write(luout,1242) lattice_unita_small(1,2), 686 > lattice_unita_small(2,2), 687 > lattice_unita_small(3,2) 688 write(luout,1243) lattice_unita_small(1,3), 689 > lattice_unita_small(2,3), 690 > lattice_unita_small(3,3) 691 write(luout,1244) lattice_unitg_small(1,1), 692 > lattice_unitg_small(2,1), 693 > lattice_unitg_small(3,1) 694 write(luout,1245) lattice_unitg_small(1,2), 695 > lattice_unitg_small(2,2), 696 > lattice_unitg_small(3,2) 697 write(luout,1246) lattice_unitg_small(1,3), 698 > lattice_unitg_small(2,3), 699 > lattice_unitg_small(3,3) 700 call lattice_small_abc_abg(a,b,c,alpha,beta,gamma) 701 write(luout,1232) a,b,c,alpha,beta,gamma 702 write(luout,1231) lattice_omega_small() 703 write(luout,1250) lattice_ecut(), 704 > ngrid_small(1),ngrid_small(2),ngrid_small(3), 705 > pack_nwave_all(2),pack_nwave(2) 706 write(luout,1251) lattice_wcut(), 707 > ngrid_small(1),ngrid_small(2),ngrid_small(3), 708 > pack_nwave_all(3),pack_nwave(3) 709 end if 710 711 call rho_symmetrizer_print(luout) 712 713 write(luout,1270) 714 if (.not.control_translation()) write(luout,1271) 715 if (.not.control_rotation()) write(luout,1272) 716 write(luout,1280) control_time_step(),control_fake_mass() 717 write(luout,1290) control_tole(),control_tolc(),control_tolr() 718 write(luout,1281) control_it_in()*control_it_out(), 719 > control_it_in(),control_it_out() 720 721 if (control_fractional()) then 722 write(6,1297) 723 if (control_fractional_smeartype().eq.-1) 724 > write(6,1298) "fixed occupation" 725 if (control_fractional_smeartype().eq.0) 726 > write(6,1298) "step function" 727 if (control_fractional_smeartype().eq.1) 728 > write(6,1298) "Fermi-Dirac" 729 if (control_fractional_smeartype().eq.2) 730 > write(6,1298) "Gaussian" 731 if (control_fractional_smeartype().eq.4) 732 > write(6,1298) "Marzari-Vanderbilt" 733 if (control_fractional_smeartype().ge.0) 734 > write(6,1299) control_fractional_kT(), 735 > control_fractional_temperature(), 736 > control_fractional_alpha() 737 end if 738 739 write(luout,1300) 740 write(luout,1305) 741 call util_flush(luout) 742 end if 743 744* |***************************| 745****************** start iterations ********************** 746* |***************************| 747 748 if (oprint) call current_second(cpu2) 749 if (oprint) call nwpw_message(14) 750 751 it_in = control_it_in() 752 it_out = control_it_out() 753 icount = 0 754 1 continue 755 icount = icount + 1 756 757 call inner_loop(ispin,neq, 758 > npack1,nfft3d,nemaxq, 759 > dcpl_mb(psi1(1)),dcpl_mb(psi2(1)),dbl_mb(dn(1)), 760 > it_in,E,deltae,deltac,deltar, 761 > dbl_mb(hml(1)),dbl_mb(lmd(1)),dbl_mb(lmd1(1)), 762 > first_iteration, 763 > dcpl_mb(psir(1)),dcpl_mb(Hpsi(1)), 764 > fractional,dbl_mb(occ1(1)),dbl_mb(occ2(1))) 765 if (fractional) E(1) = E(1) + smearcorrection 766 767 if (oprint) then 768 write(luout,1310) icount*it_in,E(1),deltae,deltac,deltar 769 call util_flush(luout) 770 end if 771 if ((deltae.gt.0.0d0).and.(icount.gt.1)) then 772 if ((icount.ge.it_out).or.(.not.fractional)) then 773 if (oprint) 774 > write(luout,*) 775 > ' *** Energy going up. iteration terminated.' 776 go to 2 777 end if 778 end if 779 deltae = dabs(deltae) 780 if ((deltae.lt.control_tole()).and. 781 > (deltac.lt.control_tolc()).and. 782 > (deltar.lt.control_tolr())) then 783 if (oprint) 784 > write(luout,*) 785 > ' *** tolerance ok. iteration terminated.' 786 go to 2 787 end if 788 789* ***** define fractional occupation **** 790 if ((icount.lt.it_out).and.fractional.and. 791 > (smeartype.ge.0)) then 792 793 call Dneall_m_diagonalize(0,dbl_mb(hml(1)), 794 > dbl_mb(eig(1)),.true.) 795c call Dneall_fmf_Multiply(0,dcpl_mb(psi2(1)),npack1, 796c > dbl_mb(hml(1)), 1.0d0, 797c > dcpl_mb(psi1(1)),0.0d0) 798 799 call psi_0define_occupation(-1.0d0,.false., 800 > ispin,ne, 801 > dbl_mb(eig(1)),dbl_mb(hml(1)), 802 > dbl_mb(occ2(1)), 803 > smeartype,smearkT, 804 > smearfermi,smearcorrection) 805c write(*,*) "definefrac smearfermi=",smearfermi 806c do i=1,ne(1) 807c write(*,*) "i,occ=",i,dbl_mb(eig(1)+i-1), 808c > dbl_mb(occ1(1)+i-1) 809c end do 810 end if 811 812 813 814 if (icount.lt.it_out) go to 1 815 if (oprint) 816 > write(luout,*) 817 > '*** arrived at the Maximum iteration. terminated.' 818 819*:::::::::::::::::::: end of iteration loop ::::::::::::::::::::::::: 820 821 2 continue 822 if (oprint) CALL nwpw_message(3) 823 if (oprint) call current_second(cpu3) 824 825* |****************************************| 826*********** produce CHECK file and diagonalize hml ***************** 827* |****************************************| 828 829* **** produce CHECK FILE **** 830 if (oprint) then 831 call util_file_name('CHECK',.true., 832 > .false., 833 > full_filename) 834 open(unit=17,file=full_filename,form='formatted') 835 end if 836 837* **** check total number of electrons **** 838 do ms =1,ispin 839 call D3dB_r_dsum(1,dbl_mb(dn(1)+(ms-1)*n2ft3d),sumall) 840 en1(ms) = sumall*lattice_omega() 841 > /dble(ngrid(1)*ngrid(2)*ngrid(3)) 842 end do 843 844 if (psp_pawexist()) then 845 if (.not.BA_push_get(mt_dbl,n2ft3d,'tmp1',tmp1(2),tmp1(1))) 846 > call errquit( 847 > 'cgsdv5: out of stack memory',0,MA_ERR) 848 849 call psp_qlm_atom(ispin,neq,dcpl_mb(psi1(1))) 850 do ms=1,ispin 851 call nwpw_compcharge_gen_dn_cmp_smooth_ms(ms,dbl_mb(tmp1(1))) 852 call Pack_c_unpack(0,dbl_mb(tmp1(1))) 853 call D3dB_cr_fft3b(1,dbl_mb(tmp1(1))) 854 call D3dB_r_Zero_Ends(1,dbl_mb(tmp1(1))) 855 call D3dB_r_dsum(1,dbl_mb(tmp1(1)),sumall) 856 en2(ms) = sumall*lattice_omega() 857 > /dble(ngrid(1)*ngrid(2)*ngrid(3)) 858 end do 859 if (.not.BA_pop_stack(tmp1(2))) 860 > call errquit( 861 > 'cgsdv5: popping stack memory',0,MA_ERR) 862 else 863 en2(1) = 0.0d0 864 en2(2) = 0.0d0 865 end if 866 en(1) = en1(1)+en2(1) 867 en(2) = en1(2)+en2(2) 868 if (omaster) then 869 write(17,1320) (en(ms),ms=1,ispin) 870 if (psp_pawexist()) then 871 write(17,1322) (en1(ms),ms=1,ispin) 872 write(17,1323) (en2(ms),ms=1,ispin) 873 end if 874 end if 875 876 877 878* **** comparison between hamiltonian an lambda matrix **** 879 if (omaster) write(17,1330) 880 do ms=1,ispin 881 do j=1,ne(ms) 882 do i=1,ne(ms) 883 w = Dneall_m_value(0,ms,i,j,dbl_mb(hml(1))) 884 virial = Dneall_m_value(0,ms,i,j,dbl_mb(lmd(1))) 885 886 if (omaster) 887 > write(17,1340) ms,i,j,w,virial,w-virial 888 end do 889 end do 890 end do 891 892 893* **** check orthonormality **** 894 if (omaster) then 895 write(17,1350) 896 end if 897 call Dneall_ffm_Multiply(0,dcpl_mb(psi1(1)), 898 > dcpl_mb(psi1(1)),npack1, 899 > dbl_mb(lmd(1))) 900 do ms=1,ispin 901 do j=1,ne(ms) 902 do i=j,ne(ms) 903 w = Dneall_m_value(0,ms,i,j,dbl_mb(lmd(1))) 904 if (omaster) write(17,1360) ms,i,j,w 905 end do 906 end do 907 end do 908 909* **** close check file **** 910 if (omaster) then 911 close(17) 912 end if 913 914* ***** do not diagonalize the hamiltonian matrix ***** 915 if (pspw_SIC().or.mparallelized) then 916 call dcopy(2*npack1*nemaxq, 917 > dcpl_mb(psi1(1)),1, 918 > dcpl_mb(psi2(1)),1) 919 920* ***** diagonalize the hamiltonian matrix but don't rotate **** 921 else if (fractional) then 922 n = ne(1) 923c call Dneall_m_HmltimesSA(0,dbl_mb(hml(1)),dbl_mb(fweight(1))) 924 call Dneall_m_diagonalize(0,dbl_mb(hml(1)), 925 > dbl_mb(eig(1)),.false.) 926c do ms=1,ispin 927c do ii=1,ne(ms) 928c dbl_mb(eig(1)+(ii-1)+(ms-1)*n) 929c > =dbl_mb(eig(1)+(ii-1)+(ms-1)*n) 930c > /dbl_mb(fweight(1)+(ii-1)+(ms-1)*n) 931c end do 932c end do 933 !*** reverse the occupation *** 934 do ms=1,ispin 935 do i=1,ne(ms) 936 dbl_mb(occ1(1)+(ms-1)*ne(1) + i-1) = 937 > dbl_mb(occ2(1)+(ms-1)*ne(1) + ne(ms)-i) 938 end do 939 end do 940 !call dcopy(ne(1)+ne(2),dbl_mb(occ2(1)),1,dbl_mb(occ1(1)),1) 941 call dcopy(2*npack1*nemaxq, 942 > dcpl_mb(psi1(1)),1, 943 > dcpl_mb(psi2(1)),1) 944 945* ***** diagonalize and rotate the hamiltonian matrix **** 946 else 947 call Dneall_m_diagonalize(0,dbl_mb(hml(1)), 948 > dbl_mb(eig(1)),.false.) 949 950* **** do not rotate if wannier localization *** 951 if (.not.pspw_HFX_localize2()) 952 > call Dneall_fmf_Multiply(0,dcpl_mb(psi1(1)),npack1, 953 > dbl_mb(hml(1)), 1.0d0, 954 > dcpl_mb(psi2(1)),0.0d0) 955 end if 956 957* |***************************| 958****************** report summary of results ********************** 959* |***************************| 960 call center_geom(cx,cy,cz) 961 call center_mass(gx,gy,gz) 962 963 if (oprint) then 964 write(luout,1300) 965 write(luout,1410) 966 write(luout,1420) 967 do I=1,ion_nion()+ion_nion_q() 968 if (I.gt.ion_nion()) then 969 write(luout,1195) I,ion_aname(I),(ion_rion(k,i),K=1,3), 970 > nwpw_cosmo_qc0(I-ion_nion()),ion_amm(i) 971 else if (ion_q_FixIon(I)) then 972 write(luout,1191) I,ion_aname(I),(ion_rion(k,i),K=1,3), 973 > ion_amass(I)/1822.89d0,ion_amm(i) 974 else if (ion_q_xyzFixIon(I)) then 975 write(6,1194) I,ion_aname(I),(ion_rion(k,i),K=1,3), 976 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 977 else 978 write(luout,1190) I,ion_aname(I),(ion_rion(k,i),K=1,3), 979 > ion_amass(I)/1822.89d0,ion_amm(i) 980 end if 981 end do 982 write(luout,1200) cx,cy,cz 983 write(luout,1210) gx,gy,gz 984 write(luout,1211) ion_nconstraints(),ion_ndof() 985 986 call pspw_charge_Print(luout) 987 988 write(luout,*) 989 write(luout,1320) en(1),en(ispin),' (real space)' 990 if (psp_pawexist()) then 991 write(luout,1322) en1(1),en1(ispin),' (real space)' 992 write(luout,1323) en2(1),en2(ispin),' (real space)' 993 end if 994 995* **** write APC potential and charges *** 996 if (pspw_V_APC_on()) call pspw_shortprint_APC(luout) 997 998 999 write(luout,1430) E(1),E(1)/ion_nion() 1000 if (pspw_qmmm_found()) then 1001 write(luout,1431) 1002 write(luout,1432) 1003 write(luout,1433) (E(1)-E(11)), 1004 > (E(1)-E(11))/ion_nion() 1005 end if 1006 1007 if (field_exist) then 1008 write(luout,1431) 1009 write(luout,1432) 1010 if (pspw_Efield_type().eq.0) then 1011 write(luout,1433) (E(1)-E(49)-E(51)), 1012 > (E(1)-E(49)-E(51))/ion_nion() 1013 else 1014 write(luout,1433) (E(1)-E(49)-E(50)-E(51)), 1015 > (E(1)-E(49)-E(50)-E(51))/ion_nion() 1016 end if 1017 end if 1018 1019 !*** print out total paw energy including core energies **** 1020 if (psp_pawexist()) 1021 > write(luout,1434) (E(1)+E(36)+E(45)), 1022 > (E(1)+E(36)+E(45))/ion_nion() 1023 1024 write(luout,1440) E(2),E(2)/n2(ispin) 1025 write(luout,1450) E(3),E(3)/n2(ispin) 1026 write(luout,1460) E(4),E(4)/n2(ispin) 1027 if (pspw_SIC()) then 1028 write(luout,1455) E(16),E(16)/n2(ispin) 1029 write(luout,1456) E(17),E(17)/n2(ispin) 1030 end if 1031 if (pspw_HFX()) then 1032 write(luout,1457) E(20),E(20)/n2(ispin) 1033 end if 1034 if (psp_U_psputerm()) then 1035 write(luout,1458) E(29),E(29)/n2(ispin) 1036 end if 1037 if (meta_found()) then 1038 write(luout,1459) E(31),E(31)/ion_nion() 1039 end if 1040 if (pspw_V_APC_on()) then 1041 write(luout,1505) E(52),E(52)/ion_nion() 1042 end if 1043 1044 1045 write(luout,1470) E(5),E(5)/ion_nion() 1046 write(luout,1480) E(6),E(6)/n2(ispin) 1047 write(luout,1490) E(7),E(7)/n2(ispin) 1048 write(luout,1495) E(8),E(8)/n2(ispin) 1049 write(luout,1496) E(9),E(9)/n2(ispin) 1050 write(luout,1497) E(10),E(10)/n2(ispin) 1051 if (pspw_SIC().and.pspw_SIC_relaxed()) then 1052 write(luout,1499) E(18),E(18)/n2(ispin) 1053 write(luout,1501) E(19),E(19)/n2(ispin) 1054 end if 1055 if (pspw_HFX().and.pspw_HFX_relaxed()) then 1056 write(luout,1502) E(21),E(21)/n2(ispin) 1057 end if 1058 if (psp_U_psputerm()) then 1059 write(luout,1503) E(30),E(30)/n2(ispin) 1060 end if 1061 if (meta_found()) then 1062 write(luout,1504) E(32),E(32)/n2(ispin) 1063 end if 1064 if (pspw_V_APC_on()) then 1065 write(luout,1506) E(53),E(53)/n2(ispin) 1066 end if 1067 1068 virial = (E(10)+E(9)+E(8)+E(7))/E(6) 1069 write(luout,1498) virial 1070 1071c **** paw energies **** 1072 if (psp_pawexist()) then 1073 write(luout,1680) E(35),E(35)/n2(ispin) 1074 write(luout,1681) E(36),E(36)/ion_nion() 1075 write(luout,1682) E(37),E(37)/n2(ispin) 1076 1077 write(luout,1683) E(38),E(38)/n2(ispin) 1078 write(luout,1684) E(39),E(39)/n2(ispin) 1079 1080 write(luout,1685) E(40),E(40)/n2(ispin) 1081 write(luout,1686) E(41),E(41)/n2(ispin) 1082 write(luout,1687) E(42),E(42)/n2(ispin) 1083 write(luout,1688) E(43),E(43)/n2(ispin) 1084 write(luout,1689) E(44),E(44)/ion_nion() 1085 write(luout,1690) E(45),E(45)/ion_nion() 1086 end if 1087 1088 if (pspw_qmmm_found()) then 1089 write(luout,1700) 1090 write(luout,1701) 1091 write(luout,1702) E(11) 1092 write(luout,1703) E(12) 1093 write(luout,1704) E(13) 1094 write(luout,1705) E(14) 1095c 1096c write(luout,1703) E(14)+E(15) 1097c write(luout,1710) E(14) 1098c write (luout,1711) E(15) 1099c 1100c write(luout,1704) E(16) 1101c write(luout,1705) E(17) 1102c write(luout,1706) E(18) 1103 end if 1104 if (ion_disp_on()) then 1105 write(luout,1720) E(33) 1106 end if 1107 1108 if (field_exist) then 1109 write(luout,1800) 1110 write(luout,1801) 1111 if (pspw_Efield_type().eq.0) then 1112 write(luout,1805) E(49)+E(51) 1113 write(luout,1806) E(49) 1114 write(luout,1807) E(50) 1115 !write(luout,1804) E(51) 1116 call dipole_Efield_print(luout) 1117 else 1118 write(luout,1805) E(49)+E(50)+E(51) 1119 write(luout,1802) E(49) 1120 write(luout,1803) E(50) 1121 write(luout,1804) E(51) 1122 end if 1123 end if 1124 1125 1126 1127 if (.not.mparallelized) then 1128* **** write out diagonal <psi|H|psi> matrix **** 1129 if (pspw_SIC()) then 1130 1131 n = ne(1) 1132 nn = n*n 1133 do ms=1,ispin 1134 if (ms.eq.1) write(luout,1331) 1135 if (ms.eq.2) write(luout,1332) 1136 !*** call Gainsville matrix output *** 1137 call output(dbl_mb(hml(1)+(ms-1)*nn), 1138 > 1,ne(ms),1,ne(ms), 1139 > n,n,1) 1140 1141 end do 1142 1143* **** write out KS eigenvalues **** 1144 else 1145 EV=27.2116d0 1146 if (control_fractional()) then 1147 if (ispin.eq.1) then 1148 write(luout,1507) smearfermi(1),smearfermi(1)*EV 1149 else 1150 write(luout,1507) smearfermi(1),smearfermi(1)*EV, 1151 > smearfermi(2),smearfermi(2)*EV 1152 end if 1153 end if 1154 write(luout,1500) 1155 NN=NE(1)-NE(2) 1156 if (fractional) then 1157 do i=1,NN 1158 write(luout,1511) dbl_mb(EIG(1)+i-1), 1159 > dbl_mb(EIG(1)+i-1)*EV, 1160 > dbl_mb(occ1(1)+i-1) 1161 end do 1162 do i=1,ne(2) 1163 write(luout,1511) dbl_mb(EIG(1)+i-1+NN), 1164 > dbl_mb(EIG(1)+i-1+NN)*EV, 1165 > dbl_mb(occ1(1)+i-1+NN), 1166 > dbl_mb(EIG(1)+i-1+n1(2)-1), 1167 > dbl_mb(EIG(1)+i-1+n1(2)-1)*EV, 1168 > dbl_mb(occ1(1)+i-1+n1(2)-1) 1169 end do 1170 else 1171 do i=1,NN 1172 write(luout,1510) dbl_mb(EIG(1)+i-1),dbl_mb(EIG(1)+i-1)*EV 1173 end do 1174 do i=1,ne(2) 1175 write(luout,1510) dbl_mb(EIG(1)+i-1+NN), 1176 > dbl_mb(EIG(1)+i-1+NN)*EV, 1177 > dbl_mb(EIG(1)+i-1+n1(2)-1), 1178 > dbl_mb(EIG(1)+i-1+n1(2)-1)*EV 1179 end do 1180 end if 1181 1182 end if 1183 end if 1184 end if 1185 1186* **** write out COSMO energies **** 1187 if (nwpw_cosmo_on()) then 1188 Egas = control_gas_energy() 1189 if (oprint) then 1190 write(luout,1730) 1191 write(luout,1731) nwpw_cosmo_screen() 1192 write(luout,1732) E(46) 1193 write(luout,1733) E(47) 1194 write(luout,1734) E(48) 1195 write(luout,1735) E(46)+E(47)+E(48), 1196 > (E(46)+E(47)+E(48))*27.2116d0*23.06d0 1197 if (dabs(Egas).gt.1.0d-6) 1198 > write(luout,1736) E(1)-Egas, 1199 > (E(1)-Egas)*27.2116d0*23.06d0 1200 call nwpw_cosmo_print(luout,Egas,E(1)) 1201 end if 1202 end if 1203 1204* **** write out extended Born solvation energies **** 1205 if (nwpw_born_on()) then 1206 Egas = control_gas_energy() 1207 if (oprint) then 1208 write(luout,1740) 1209 write(luout,1741) nwpw_born_screen() 1210 write(luout,1745) E(52),E(52)*27.2116d0*23.06d0 1211 if (dabs(Egas).gt.1.0d-6) 1212 > write(luout,1746) E(1)-Egas, 1213 > (E(1)-Egas)*27.2116d0*23.06d0 1214 call nwpw_born_print(luout,Egas,E(1)) 1215 end if 1216 end if 1217 1218 1219 if (mparallelized) then 1220 do ms=1,ispin 1221 if ((ms.eq.1).and.(oprint)) write(luout,1333) 1222 if ((ms.eq.2).and.(oprint)) write(luout,1334) 1223 do i=1,ne(ms) 1224 w = Dneall_m_value(0,ms,i,i,dbl_mb(hml(1))) 1225 if (oprint) write(luout,1510) w,w*EV 1226 end do 1227 end do 1228 end if 1229 1230 if (oprint) then 1231* ***** extra energy output for QA test **** 1232 write(luout,1600) E(1) 1233 end if 1234 1235 1236 1237* |***************************| 1238****************** Prologue ********************** 1239* |***************************| 1240 1241* **** calculate spin contamination **** 1242 call Calculate_psi_spin2(ispin,ne,npack1,dcpl_mb(psi2(1)), 1243 > fractional,dbl_mb(occ2(1)),w) 1244 1245* **** calculate the Dipole *** 1246 call Calculate_Dipole(ispin,ne,n2ft3d,dbl_mb(dn(1)),dipole) 1247 1248* **** perfom Lubin and Mulliken analysis *** 1249 if (control_Mulliken()) then 1250 1251* **** Lubin Water Analysis *** 1252 call pspw_Lubin_water_analysis(rtdb,ispin,ne,n2ft3d, 1253 > dbl_mb(dn(1))) 1254 1255* **** Analysis *** 1256 call pspw_analysis(0,rtdb,ispin,ne,dcpl_mb(psi2(1)), 1257 > dbl_mb(eig(1))) 1258 1259* **** generate APC ***** 1260 call pspw_dngen_APC(ispin,ne,dbl_mb(dn(1)),.false.) 1261 call pspw_print_APC(luout) 1262 1263 end if 1264 1265* ***** write psi2 wavefunctions **** 1266 call psi_write(ispin,ne,dcpl_mb(psi2(1)), 1267 > smearoccupation,dbl_mb(occ2(1))) 1268 1269* **** write geometry to rtdb **** 1270 call pspw_charge_write(rtdb) 1271 call ion_write(rtdb) 1272 1273 1274* **** deallocate heap memory **** 1275 if (control_version().eq.3) call ewald_end() 1276 call strfac_end() 1277 if (control_version().eq.3) call coulomb_end() 1278 if (control_version().eq.4) call coulomb2_end() 1279 call ke_end() 1280 call mask_end() 1281 call Pack_end() 1282 call G_end() 1283 call psp_U_end() 1284 call vdw_DF_end() 1285 call nwpw_meta_gga_end() 1286 call rho_symmetrizer_end() 1287 call pspw_end_SIC() 1288 call pspw_end_HFX() 1289 call pspw_end_APC() 1290 call pspw_qmmm_end() 1291 call pspw_charge_end() 1292 call pspw_Efield_end() 1293 call meta_finalize(rtdb) 1294 if (psp_pawexist()) call nwpw_gintegrals_end() 1295 1296 call ion_end() 1297 call psp_end() 1298 call ion_end_FixIon() 1299 1300c call frac_occ_end() 1301 1302 value = BA_free_heap(psir(2)) 1303 value = BA_free_heap(Hpsi(2)) 1304 value = BA_free_heap(dn(2)) 1305 value = BA_free_heap(eig(2)) 1306 value = Dneall_m_free(hml) 1307 value = Dneall_m_free(lmd) 1308 value = Dneall_m_free(lmd1) 1309 value = BA_free_heap(psi1(2)) 1310 value = BA_free_heap(psi2(2)) 1311 if (fractional) then 1312 value = BA_free_heap(occ1(2)) 1313 value = BA_free_heap(occ2(2)) 1314 end if 1315 call D3dB_pfft_end() 1316 call D3dB_end(1) 1317 if (control_version().eq.4) call D3dB_end(2) 1318 if (control_has_ngrid_small()) call D3dB_end(3) 1319 call Dne_end() 1320 call psi_data_end() 1321 1322* |***************************| 1323****************** report consumed cputime ********************** 1324* |***************************| 1325 if (oprint) then 1326 CALL current_second(cpu4) 1327 1328 T1=CPU2-CPU1 1329 T2=CPU3-CPU2 1330 T3=CPU4-CPU3 1331 T4=CPU4-CPU1 1332 AV=T2/dble(icount*it_in) 1333 write(luout,*) 1334 write(luout,*) '-----------------' 1335 write(luout,*) 'cputime in seconds' 1336 write(luout,*) 'prologue : ',T1 1337 write(luout,*) 'main loop : ',T2 1338 write(luout,*) 'epilogue : ',T3 1339 write(luout,*) 'total : ',T4 1340 write(luout,*) 'cputime/step: ',AV 1341 write(luout,*) 1342 call nwpw_timing_print_final(.true.,(icount*it_in)) 1343 CALL nwpw_MESSAGE(4) 1344 end if 1345 1346 1347 call Parallel2d_Finalize() 1348 call Parallel_Finalize() 1349 cpsdv5 = value 1350 return 1351 1352 1353*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 1354 1000 FORMAT(10X,'****************************************************') 1355 1010 FORMAT(10X,'* *') 1356 1020 FORMAT(10X,'* Car-Parrinello microcluster calculation *') 1357 1030 FORMAT(10X,'* [ steepest descent minimization ] *') 1358 1035 FORMAT(10x,'* [ NorthWest Chemistry implementation ] *') 1359 1040 FORMAT(10X,'* version #5.00 06/01/99 *') 1360 1041 FORMAT(10X,'* This code was developed by Eric J. Bylaska, *') 1361 1042 FORMAT(10X,'* and was based upon algorithms and code *') 1362 1043 FORMAT(10X,'* developed by the group of Prof. John H. Weare *') 1363 1100 FORMAT(//) 1364 1110 FORMAT(10X,'================ PSPW input data ===================') 1365 1111 FORMAT(/' number of processors used:',I10) 1366 1112 FORMAT( ' parallel mapping : 1d-slab') 1367 1113 FORMAT( ' parallel mapping : 2d-hilbert') 1368 1114 FORMAT( ' parallel mapping : balanced') 1369 1115 FORMAT(/' options:') 1370 1116 FORMAT( ' parallel mapping : not balanced') 1371 1117 FORMAT( ' processor grid :',I4,' x',I4) 1372 1118 FORMAT( ' parallel mapping : 2d-hcurve') 1373 1119 FORMAT( ' parallel io : on') 1374 1120 FORMAT(5X,' ionic motion = ',A) 1375 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 1376 1122 FORMAT( ' parallel io : off') 1377 1123 FORMAT( ' number of threads :',I10) 1378 1130 FORMAT(5X,' electron spin = ',A) 1379 1131 FORMAT(5X,' exchange-correlation = ',A) 1380 1132 FORMAT(5X,' using fractional occupation') 1381 1135 FORMAT(/' The masses of QM H atoms converted to 2.0 amu. ', 1382 > /' To turn off this default', 1383 > /' nwpw', 1384 > /' makehmass2 off', 1385 > /' end') 1386 1140 FORMAT(/' elements involved in the cluster:') 1387 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I1) 1388 1151 FORMAT(5X,' cutoff =',4F8.3) 1389 1152 FORMAT(12X,' highest angular component : ',i2) 1390 1153 FORMAT(12X,' local potential used : ',i2) 1391 1154 FORMAT(12X,' number of non-local projections: ',i2) 1392 1155 FORMAT(12X,' semicore corrections included : ', 1393 > F6.3,' (radius) ',F6.3,' (charge)') 1394 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 1395 1159 FORMAT(/' total charge=',F8.3) 1396 1160 FORMAT(/' atomic composition:') 1397 1170 FORMAT(7(5X,A2,':',I5)) 1398 1180 FORMAT(/' initial position of ions:') 1399 1190 FORMAT(5X, I4, A5, ' (',3F11.5,' ) - atomic mass= ',F7.3,' ',A) 1400 1191 FORMAT(5X, I4, A5, ' (',3F11.5, 1401 > ' ) - atomic mass= ',F7.3,' - fixed ',A) 1402 1193 FORMAT(5X, I4, A5, ' (',3F11.5, 1403 > ' ) - atomic mass= ',F7.3,' - z fixed') 1404 1194 FORMAT(5X, I4, A5, ' (',3F11.5, 1405 > ' ) - atomic mass= ',F7.3,A) 1406 1195 FORMAT(5X, I4, A5, ' (',3F11.5,' ) - charge= ',E13.6,' ',A) 1407 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 1408 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 1409 1211 FORMAT(5X,' number of constraints = ', I6,' ( DOF = ',I6,' )' ) 1410 1219 FORMAT(/' number of electrons: spin up=',F6.2, 16x, 1411 > ' down=',F6.2,A) 1412 1220 FORMAT(/' number of electrons: spin up=',I6, 1413 > ' (',I4,' per task)', 1414 > ' down=',I6, 1415 > ' (',I4,' per task)', 1416 > A) 1417 1221 FORMAT( ' number of orbitals : spin up=',I6, 1418 > ' (',I4,' per task)', 1419 > ' down=',I6, 1420 > ' (',I4,' per task)', 1421 > A) 1422 1229 FORMAT(/' small supercell:') 1423 1230 FORMAT(/' supercell:') 1424 1231 FORMAT(5x,' volume : ',F12.1) 1425 1232 FORMAT(5x,' lattice: a= ',f8.3,' b= ',f8.3,' c= ',f8.3, 1426 > /5x,' alpha=',f8.3,' beta=',f8.3,' gamma=',f8.3) 1427 1233 FORMAT(5x,' cell_name: ',A) 1428 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 1429 1242 FORMAT(5x,' a2=<',3f8.3,' >') 1430 1243 FORMAT(5x,' a3=<',3f8.3,' >') 1431 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >') 1432 1245 FORMAT(5x,' b2=<',3f8.3,' >') 1433 1246 FORMAT(5x,' b3=<',3f8.3,' >') 1434 1435 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1436 & '( ',I8,' waves ',I8,' per task)') 1437 1251 FORMAT(5X,' wavefnc cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1438 & '( ',I8,' waves ',I8,' per task)') 1439 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,' and',I3) 1440 1261 FORMAT(5X,' madelung=',f14.8) 1441 1270 FORMAT(/' technical parameters:') 1442 1271 FORMAT(5x, ' translation constrained') 1443 1272 FORMAT(5x, ' rotation constrained') 1444 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1) 1445 1281 FORMAT(5X, ' maximum iterations =',I10, 1446 > ' ( ',I4,' inner ',I6,' outer )') 1447 1290 FORMAT(5X, ' tolerance=',E9.3,' (energy)',E12.3, 1448 & ' (electron)',E12.3,' (ion)') 1449 1297 FORMAT(/' fractional smearing parameters:') 1450 1298 FORMAT(5X, ' smearing algorithm = ',A) 1451 1299 FORMAT(5X, ' smearing parameter = ',E9.3,' (',F7.1,' K)'/, 1452 > 5X, ' mixing parameter =',F7.4) 1453 1300 FORMAT(//) 1454 1305 FORMAT(10X, 1455 > '==================== iteration =========================') 1456 1310 FORMAT(I8,E20.10,3E15.5) 1457 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 1458 1321 FORMAT(' total charge of system:',F11.5,A) 1459 1322 FORMAT(' plane-wave part: ',F11.5,' ',F11.5,A) 1460 1323 FORMAT(' augmented part: ',F11.5,' ',F11.5,A) 1461 1462 1330 FORMAT(/' comparison between hamiltonian and lambda matrix') 1463 1331 FORMAT(/' Elements of Hamiltonian matrix (up/restricted)') 1464 1332 FORMAT(/' Elements of Hamiltonian matrix (down)') 1465 1333 FORMAT(/' Diagonal Elements of Hamiltonian matrix(up/restricted)') 1466 1334 FORMAT(/' Diagonal Elements of Hamiltonian matrix(down)') 1467 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7) 1468 1341 FORMAT(I3,2I3,' H=',E16.6) 1469 1350 FORMAT(/' orthonormality') 1470 1360 FORMAT(I3,2I3,E18.7) 1471 1370 FORMAT(I3) 1472 1380 FORMAT(' ''',a,'''',I4) 1473 1390 FORMAT(I3) 1474 1400 FORMAT(I3,3E18.8/3X,3E18.8) 1475 1410 FORMAT(10X,'============= summary of results =================') 1476 1420 FORMAT( ' final position of ions:') 1477 1430 FORMAT(//' total energy :',E19.10,' (',E15.5,'/ion)') 1478 1431 FORMAT(/' QM Energies') 1479 1432 FORMAT( '------------') 1480 1433 FORMAT( ' total QM energy :',E19.10,' (',E15.5,'/ion)') 1481 1434 FORMAT(//' total paw energy :',E19.10,' (',E15.5,'/ion)') 1482 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)') 1483 1450 FORMAT( ' hartree energy :',E19.10,' (',E15.5,'/electron)') 1484 1455 FORMAT( ' SIC-hartree energy :',E19.10,' (',E15.5,'/electron)') 1485 1456 FORMAT( ' SIC-exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1486 1457 FORMAT( ' HF exchange energy :',E19.10,' (',E15.5,'/electron)') 1487 1458 FORMAT( ' DFT+U energy :',E19.10,' (',E15.5,'/electron)') 1488 1459 FORMAT( ' Metadynamics energy :',E19.10,' (',E15.5,'/ion)') 1489 1460 FORMAT( ' exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1490 1470 FORMAT( ' ion-ion energy :',E19.10,' (',E15.5,'/ion)') 1491 1480 FORMAT(/' K.S. kinetic energy :',E19.10,' (',E15.5,'/electron)') 1492 1490 FORMAT( ' K.S. V_l energy :',E19.10,' (',E15.5,'/electron)') 1493 1495 FORMAT( ' K.S. V_nl energy :',E19.10,' (',E15.5,'/electron)') 1494 1496 FORMAT( ' K.S. V_Hart energy :',E19.10,' (',E15.5,'/electron)') 1495 1497 FORMAT( ' K.S. V_xc energy :',E19.10,' (',E15.5,'/electron)') 1496 1498 FORMAT( ' Virial Coefficient :',E19.10) 1497 1499 FORMAT( ' K.S. SIC-hartree energy :',E19.10, 1498 > ' (',E15.5,'/electron)') 1499 1500 FORMAT(/' orbital energies:') 1500 1501 FORMAT( ' K.S. SIC-exc-corr energy :',E19.10, 1501 > ' (',E15.5,'/electron)') 1502 1502 FORMAT( ' K.S. HFX energy :',E19.10, 1503 > ' (',E15.5,'/electron)') 1504 1503 FORMAT( ' K.S. DFT+U energy :',E19.10,' (',E15.5,'/electron)') 1505 1504 FORMAT( ' K.S. Metadynamics energy :',E19.10,' (',E15.5,'/ion)') 1506 1505 FORMAT( ' APC energy :',E19.10,' (',E15.5,'/ion)') 1507 1506 FORMAT( ' K.S. V_APC energy :',E19.10,' (',E15.5,'/ion)') 1508 1507 FORMAT(/' Fermi energy =',2(E18.7,' (',F8.3,'eV)')) 1509 1510 1511 1510 FORMAT(2(E18.7,' (',F8.3,'eV)')) 1512 1511 FORMAT(2(E18.7,' (',F8.3,'eV) occ=',F5.3)) 1513 1600 FORMAT(/' Total PSPW energy :',E19.10) 1514 1515 1680 FORMAT(/' kinetic (loc. basis):',E19.10,' (',E15.5,'/electron)') 1516 1681 FORMAT( ' kinetic (core):',E19.10,' (',E15.5,'/ion)') 1517 1682 FORMAT( ' valence-core :',E19.10,' (',E15.5,'/electron)') 1518 1683 FORMAT( ' ncmp*V_local :',E19.10,' (',E15.5,'/electron)') 1519 1684 FORMAT( ' V_lpaw (loc. basis):',E19.10,' (',E15.5,'/electron)') 1520 1685 FORMAT( ' coulomb (loc. basis):',E19.10,' (',E15.5,'/electron)') 1521 1686 FORMAT( ' coulomb (ncmp-ncmp):',E19.10,' (',E15.5,'/electron)') 1522 1687 FORMAT( ' coulomb (ncmp-pw):',E19.10,' (',E15.5,'/electron)') 1523 1688 FORMAT( ' exc-cor (loc. basis):',E19.10,' (',E15.5,'/electron)') 1524 1689 FORMAT( ' pxc-cor (loc. basis):',E19.10,' (',E15.5,'/electron)') 1525 1690 FORMAT( ' ion-core :',E19.10,' (',E15.5,'/ion)') 1526 1691 FORMAT( ' V_local PAW residual:',E19.10,' (',E15.5,'/electron)') 1527 1528 1700 FORMAT(/' QM/MM-pol-vib/CAV Energies') 1529 1701 FORMAT( ' --------------------------') 1530 1702 FORMAT( ' LJ energy :',E19.10) 1531 1703 FORMAT( ' Residual Coulomb energy:',E19.10) 1532 1704 FORMAT( ' MM Vibration energy :',E19.10) 1533 1705 FORMAT( ' QM/MM coupling energy :',E19.10) 1534 1706 FORMAT( ' (QM+MM)/Cavity energy :',E19.10) 1535 1536 1720 FORMAT(/' Dispersion energy :',E19.10) 1537 1538 1730 FORMAT(/' COSMO energies:') 1539 1731 FORMAT(5x,' screen=(epsilon-1)/(epsilon+1/2):',F11.6) 1540 1732 FORMAT(5x,' screen*Qelc*B*q energy :',E19.10) 1541 1733 FORMAT(5x,' screen*Qion*B*q energy :',E19.10) 1542 1734 FORMAT(5x,' screen*0.5*q*A*q energy :',E19.10) 1543 1735 FORMAT(5x,' solvation energy (w/o QM polarization) :',E19.10, 1544 > ' (',F8.3,' kcal/mol)') 1545 1736 FORMAT(5x,' solvation energy (w/ QM polarization) :',E19.10, 1546 > ' (',F8.3,' kcal/mol)') 1547 1548 1740 FORMAT(/' extended Born solvation energies:') 1549 1741 FORMAT(5x,' screen=(epsilon-1)/(epsilon):',F11.6) 1550 1745 FORMAT(5x,' solvation energy (w/o QM polarization) :',E19.10, 1551 > ' (',F8.3,' kcal/mol)') 1552 1746 FORMAT(5x,' solvation energy (w/ QM polarization) :',E19.10, 1553 > ' (',F8.3,' kcal/mol)') 1554 1555 1556 1800 FORMAT(/' Charge+Electric Field Energies') 1557 1801 FORMAT( ' ------------------------------') 1558 1802 FORMAT( ' - Charge+Electric Field/Electron :',E19.10) 1559 1803 FORMAT( ' - Charge+Electric Field/Ion :',E19.10) 1560 1804 FORMAT( ' - Charge+Electric Field/Charge Field:',E19.10) 1561 1805 FORMAT( ' Charge+Electric Field Energy :',E19.10) 1562 1806 FORMAT( ' - Electric Field/Resta_Dipole :',E19.10) 1563 1807 FORMAT( ' - K.S. Electric Field/Resta_Dipole :',E19.10) 1564 1565 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 1566 1567 9000 if (omaster) write(luout,9010) ierr 1568 call Parallel2d_Finalize() 1569 call Parallel_Finalize() 1570 1571 cpsdv5 = value 1572 return 1573 END 1574