1* 2* $Id$ 3* 4 5*********************************************************************** 6* cpmdv5-mpi (MPI code) * 7* * 8* This is a developing cpsdv5 parallel code for NWChem * 9* + mpi message passing library used * 10* + ngp is used instead of nfft in this proceudure * 11* + error checking is based on aimd.h parameters * 12* then control file * 13* + my own slap-decomposed parallel 3d-FFT(real->complex) used * 14* * 15* * 16*********************************************************************** 17 18 logical function cpmdv5(rtdb) 19CDEC$ OPTIMIZE:0 20 implicit none 21 integer rtdb 22 23#include "global.fh" 24#include "bafdecls.fh" 25#include "btdb.fh" 26#include "stdio.fh" 27#include "errquit.fh" 28cccc#include "frac_occ.fh" 29 30 logical value 31 32 real*8 kb 33 parameter (kb=3.16679d-6) 34 real*8 autoatm 35 parameter (autoatm =290.360032539d6) 36 37 38 39* **** parallel variables **** 40 integer taskid,np,np_i,np_j 41 integer MASTER 42 parameter(MASTER=0) 43 44* **** timing variables **** 45 real*8 cpu1,cpu2,cpu3,cpu4 46 real*8 t1,t2,t3,t4,av 47 48* **** lattice variables **** 49 integer ngrid(3),nwave,nfft3d,n2ft3d,ngrid_small(3) 50 integer npack1 51 real*8 a,b,c,alpha,beta,gamma 52 53* **** electronic variables **** 54 real*8 icharge 55 integer ispin 56 integer ne(2),n1(2),n2(2),nemax,neall,neq(2),nemaxq 57 real*8 en(2),en1(2),en2(2) 58 real*8 dipole(3) 59 60* complex*16 psi1(nfft3d,nemax) 61* complex*16 psi2(nfft3d,nemax) 62* real*8 dn(n2ft3d,2) 63* complex*16 Hpsi(nfft3d,nemax) 64* complex*16 psir(nfft3d,nemax) 65 integer psi0(2),psi1(2),psi2(2) 66 integer occ0(2),occ1(2),occ2(2) 67 integer dn(2) 68 integer Hpsi(2),psir(2) 69 70 71* ***** energy variables **** 72 real*8 E(60),eke,eave,evar,cv,have,hvar,qave,qvar,Egas 73 74* real*8 eig(2*nemax) 75* real*8 hml(2*nemax*nemax) 76* real*8 lmd(2*nemax*nemax) 77 integer eig(2),hml(2),lmd(2),lmd1(2) 78 79* **** psi smearing block **** 80 logical fractional 81 integer smearoccupation,smeartype 82 real*8 smearfermi(2),smearcorrection,smearkT 83 84* **** error variables **** 85 integer ierr 86 87* **** local variables **** 88 logical verlet,mulliken,SA,found,calc_pressure,found_bak 89 logical field_exist 90 integer ms 91 real*8 gx,gy,gz,cx,cy,cz 92 real*8 vgx,vgy,vgz,vcx,vcy,vcz 93 real*8 ekg,eki0,eki1,sum 94 real*8 eke0,eke1 95 real*8 EV,pi,dt 96 real*8 emotion_time_shift 97 integer i,j,k,ia,n,nn 98 integer ii,jj,index,indx 99 integer icount,it_in,it_out,icount_shift 100 real*8 w,sumall,pressure,stress(3,3),p1,p2 101 real*8 Te_init,Tr_init,Te_new,Tr_new,sa_decay(2),sa_alpha(2) 102 double precision tollz 103 parameter(tollz=1d-16) 104 integer mapping,mapping1d 105 character*50 filename 106 character*255 full_filename,full_bak 107 integer tmp1(2) 108 109 110 111* **** external functions **** 112 real*8 psp_zv,psp_rc,ewald_rcut,ion_amass 113 real*8 ewald_mandelung,lattice_omega_small 114 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 115 real*8 lattice_unitg,lattice_unitg_small,lattice_unita_small 116 integer ewald_ncut,ewald_nshl3d 117 integer psp_lmmax,psp_lmax,psp_locp,ion_nkatm 118 character*4 ion_atom,ion_aname 119 external psp_zv,psp_rc,ewald_rcut,ion_amass 120 external ewald_mandelung,lattice_omega_small 121 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 122 external lattice_unitg,lattice_unitg_small,lattice_unita_small 123 external ewald_ncut,ewald_nshl3d 124 external psp_lmmax,psp_lmax,psp_locp,ion_nkatm 125 external ion_atom,ion_aname 126 127 real*8 control_rti,control_rte,ion_rion 128 real*8 ion_vion,ion_com_ke,ion_ke 129 real*8 ion_Temperature,ion_com_Temperature 130 external control_rti,control_rte,ion_rion 131 external ion_vion,ion_com_ke,ion_ke 132 external ion_Temperature,ion_com_Temperature 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,ion_nion,ion_natm,ion_katm 145 external pack_nwave_all 146 external control_it_in,control_it_out,control_gga,control_version 147 external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm 148 149 character*12 control_boundry 150 external control_boundry 151 152 logical psp_semicore,pspw_qmmm_found,pspw_charge_found 153 external psp_semicore,pspw_qmmm_found,pspw_charge_found 154 real*8 psp_rcore,psp_ncore,psp_rlocal 155 external psp_rcore,psp_ncore,psp_rlocal 156 logical pspw_Efield_found 157 external pspw_Efield_found 158 159 logical control_Nose,control_Mulliken,Nose_restart 160 external control_Nose,control_Mulliken,Nose_restart 161 162 integer Nose_Mchain,Nose_Nchain 163 external Nose_Mchain,Nose_Nchain 164 165 real*8 control_Nose_Te,Nose_Qe,Nose_Pe,Nose_Ee0 166 external control_Nose_Te,Nose_Qe,Nose_Pe,Nose_Ee0 167 168 real*8 control_Nose_Tr,Nose_Qr,Nose_Pr,Nose_Er0 169 external control_Nose_Tr,Nose_Qr,Nose_Pr,Nose_Er0 170 logical v_psi_filefind 171 external v_psi_filefind 172 real*8 nwpw_timing 173 external nwpw_timing 174 175 logical control_out_of_time,control_new_vpsi 176 external control_out_of_time,control_new_vpsi 177 178 logical pspw_HFX_localize2 179 logical control_SA,control_Fei,pspw_SIC,pspw_HFX,control_pressure 180 real*8 control_SA_decay 181 external pspw_HFX_localize2 182 external control_SA,control_Fei,pspw_SIC,pspw_HFX,control_pressure 183 external control_SA_decay 184 185 integer control_np_orbital,control_mapping,control_mapping1d 186 external control_np_orbital,control_mapping,control_mapping1d 187 188 189 logical control_translation,control_rotation,control_balance 190 external control_translation,control_rotation,control_balance 191 192 logical Dneall_m_allocate,Dneall_m_free,control_parallel_io 193 external Dneall_m_allocate,Dneall_m_free,control_parallel_io 194 195 real*8 Dneall_m_value,pspw_qmmm_lambda 196 external Dneall_m_value,pspw_qmmm_lambda 197 198 logical meta_found,tamd_found,psp_U_psputerm,ion_disp_on 199 external meta_found,tamd_found,psp_U_psputerm,ion_disp_on 200 201 integer ion_nconstraints,ion_ndof,control_ngrid_small 202 external ion_nconstraints,ion_ndof,control_ngrid_small 203 204 logical psp_pawexist,ion_makehmass2,control_has_ngrid_small 205 external psp_pawexist,ion_makehmass2,control_has_ngrid_small 206 207 logical nwpw_born_on 208 external nwpw_born_on 209 real*8 nwpw_born_screen 210 external nwpw_born_screen 211 logical pspw_V_APC_on 212 external pspw_V_APC_on 213 real*8 control_gas_energy 214 external control_gas_energy 215 216 character*50 control_cell_name 217 external control_cell_name 218 integer Parallel_maxthreads 219 external Parallel_maxthreads 220 221 222* |************| 223*****************************| PROLOGUE |**************************** 224* |************| 225 226 value = .true. 227 pi = 4.0d0*datan(1.0d0) 228 229 call nwpw_timing_init() 230 call dcopy(60,0.0d0,0,E,1) 231 232 233* **** get parallel variables **** 234 call Parallel_Init() 235 call Parallel_np(np) 236 call Parallel_taskid(taskid) 237 if (taskid.eq.MASTER) call current_second(cpu1) 238 239* ***** print out header **** 240 if (taskid.eq.MASTER) then 241 write(luout,1000) 242 write(luout,1010) 243 write(luout,1020) 244 write(luout,1010) 245 write(luout,1030) 246 write(luout,1031) 247 write(luout,1010) 248 write(luout,1035) 249 write(luout,1010) 250 write(luout,1040) 251 write(luout,1010) 252 write(luout,1041) 253 write(luout,1042) 254 write(luout,1043) 255 write(luout,1010) 256 write(luout,1000) 257 call nwpw_message(1) 258 write(luout,1110) 259 call util_flush(luout) 260 end if 261 call ga_sync() 262 263 value = control_read(2,rtdb) 264 call Parallel2d_Init(control_np_orbital()) 265 call Parallel2d_np_i(np_i) 266 call Parallel2d_np_j(np_j) 267 268 ngrid(1) = control_ngrid(1) 269 ngrid(2) = control_ngrid(2) 270 ngrid(3) = control_ngrid(3) 271 nwave = 0 272 mapping = control_mapping() 273 274* **** initialize psi_data **** 275 call psi_data_init(100) 276 277 278* **** initialize D3dB data structure **** 279 call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 280 call D3dB_nfft3d(1,nfft3d) 281 n2ft3d = 2*nfft3d 282 if (control_version().eq.4) 283 > call D3dB_Init(2,2*ngrid(1),2*ngrid(2),2*ngrid(3),mapping) 284 285 if (control_has_ngrid_small()) then 286 ngrid_small(1) = control_ngrid_small(1) 287 ngrid_small(2) = control_ngrid_small(2) 288 ngrid_small(3) = control_ngrid_small(3) 289 call D3dB_Init(3,ngrid_small(1),ngrid_small(2),ngrid_small(3), 290 > mapping) 291 end if 292 293 294* **** initialize lattice data structure **** 295 call lattice_init() 296 call G_init() 297 call mask_init() 298 call Pack_init() 299 call Pack_npack(1,npack1) 300 301 call D3dB_pfft_init() 302 call ga_sync() 303 304* ***** allocate psi2, psi1, and psi0 wavefunctions **** 305 call psi_get_ne_occupation(ispin,ne,smearoccupation) 306 if (smearoccupation.gt.0) then 307 fractional = .true. 308 else 309 fractional = .false. 310 end if 311 mapping1d = control_mapping1d() 312 call Dne_init(ispin,ne,mapping1d) 313 call Dneall_neq(neq) 314 nemaxq = neq(1)+neq(2) 315 316 value = BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 317 > 'psi2',psi2(2),psi2(1)) 318 value = value.and. 319 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 320 > 'psi1',psi1(2),psi1(1)) 321 value = value.and. 322 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 323 > 'psi0',psi0(2),psi0(1)) 324 if (fractional) then 325 value = value.and. 326 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ0',occ0(2),occ0(1)) 327 value = value.and. 328 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ1',occ1(2),occ1(1)) 329 value = value.and. 330 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ2',occ2(2),occ2(1)) 331 end if 332 if (.not.value) call errquit('out of heap memory',0, MA_ERR) 333 334 335 336* ***** read psi2 wavefunctions **** 337 call psi_read(ispin,ne,dcpl_mb(psi2(1)), 338 > smearoccupation,dbl_mb(occ2(1))) 339 340 341* **** move wavefunction velocities **** 342 if (control_new_vpsi()) then 343 call v_psi_delete() 344 end if 345 346 347* **** generate initial wavefunction velocities if it does not exist **** 348 if (.not.v_psi_filefind()) then 349 call v_psi_new(ispin,ne) 350 end if 351 352 353* ***** read psi0 wavefunctions **** 354 call dcopy(2*npack1*(neq(1)+neq(2)),0.0d0,0,dcpl_mb(psi1(1)),1) 355 call v_psi_read(ispin,ne,dcpl_mb(psi1(1))) 356 n1(1) = 1 357 n2(1) = ne(1) 358 n1(2) = ne(1)+1 359 n2(2) = ne(1)+ne(2) 360 nemax = ne(1)+ne(2) 361 362* **** allocate other variables ***** 363 value = BA_alloc_get(mt_dbl,(2*nemax),'eig',eig(2),eig(1)) 364 value = value.and.Dneall_m_allocate(0,hml) 365 value = value.and.Dneall_m_allocate(0,lmd) 366 value = value.and.Dneall_m_allocate(0,lmd1) 367 368 value = value.and. 369 > BA_alloc_get(mt_dbl,(4*nfft3d), 370 > 'dn',dn(2),dn(1)) 371 value = value.and. 372 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 373 > 'Hpsi',Hpsi(2),Hpsi(1)) 374 value = value.and. 375 > BA_alloc_get(mt_dcpl,nfft3d*(neq(1)+neq(2)), 376 > 'psir',psir(2),psir(1)) 377 if (.not.value) call errquit('out of heap memory',0, MA_ERR) 378 379 380 381* **** read ions **** 382 value = ion_init(rtdb) 383 384 385* **** initialize FixIon constraint **** 386 call ion_init_FixIon(rtdb) 387 388 389* **** allocate psp data structure and read in psedupotentials into it **** 390 call psp_init() 391 call psp_readall() 392 if (psp_semicore(0)) call semicore_check() 393 394 395* **** initialize G,mask,ke,and coulomb data structures **** 396 call ke_init() 397 if (control_version().eq.3) call coulomb_init() 398 if (control_version().eq.4) call coulomb2_init() 399 call strfac_init() 400 if (control_version().eq.3) call ewald_init() 401 402* **** initialize QM/MM **** 403 call pspw_init_APC(rtdb) 404 call pspw_qmmm_init(rtdb) 405 call pspw_charge_init(rtdb) 406 call pspw_Efield_init(rtdb,ispin,ne) 407 field_exist = pspw_charge_found().or.pspw_Efield_found() 408 409* ****************************** 410* **** scaling psi velocity **** 411* ****************************** 412 call dcopy(2*(neq(1)+neq(2))*npack1,dcpl_mb(psi1(1)),1, 413 > dcpl_mb(psi0(1)),1) 414 call dscal(2*(neq(1)+neq(2))*npack1,control_rte(), 415 > dcpl_mb(psi1(1)),1) 416 eke0 = 0.0d0 417 eke1 = 0.0d0 418 do i=1,(neq(1)+neq(2)) 419 call Pack_cc_dot(1,dcpl_mb(psi0(1)+(i-1)*npack1), 420 > dcpl_mb(psi0(1)+(i-1)*npack1), 421 > sum) 422 eke0 = eke0 + sum 423 call Pack_cc_dot(1,dcpl_mb(psi1(1)+(i-1)*npack1), 424 > dcpl_mb(psi1(1)+(i-1)*npack1), 425 > sum) 426 eke1 = eke1 + sum 427 end do 428 429 call D1dB_SumAll(eke0) 430 call D1dB_SumAll(eke1) 431 eke0 = control_fake_mass()*eke0 432 eke1 = control_fake_mass()*eke1 433 call ion_init_ke(ekg,eki0,eki1) 434 435 436* **** Initialize thermostats **** 437 if (control_Nose()) then 438 call ke_ave(ispin,neq,dcpl_mb(psi2(1)),w, 439 > fractional,dbl_mb(occ2(1))) 440 call Nose_Init((ne(1)+ne(2)),w) 441 end if 442 443 444* **** Initialize simulated annealing **** 445 SA = .false. 446 Te_init = 0.0d0 447 Tr_init = 0.0d0 448 sa_alpha(1) = 1.0d0 449 sa_alpha(2) = 1.0d0 450 if (control_SA()) then 451 if (control_Nose()) then 452 SA = .true. 453 sa_decay(1) = control_SA_decay(1) 454 sa_decay(2) = control_SA_decay(2) 455 Te_init = control_Nose_Te() 456 Tr_init = control_Nose_Tr() 457 else 458 dt = control_time_step() 459 SA = .false. 460 sa_decay(1) = control_SA_decay(1) 461 sa_decay(2) = control_SA_decay(2) 462 sa_alpha(1) = dexp( -(dt/control_SA_decay(1)) ) 463 sa_alpha(2) = dexp( -(dt/control_SA_decay(2)) ) 464 end if 465 end if 466 467 468* **** initialize two-electron Gaussian integrals **** 469* **** initialize paw ncmp*Vloc **** 470 if (psp_pawexist()) then 471 call nwpw_gintegrals_init() 472 call nwpw_gintegrals_set(.true.) 473 call psp_dE_ncmp_vloc_Qlm(ispin,.false.,dipole) 474 end if 475 476* **** initialize metadynamics and tamd **** 477 call meta_initialize(rtdb) 478 call tamd_initialize(rtdb) 479 480* **** initialize dplot **** 481 call dplot_iteration_init() 482 483c* **** initialize frac_occ data structure **** 484c call frac_occ_init(rtdb,ispin,ne) 485 486* **** initialize SIC and HFX **** 487 call pspw_init_SIC(rtdb,ne) 488 call pspw_init_HFX(rtdb,ispin,ne) 489 490 491* **** initialize DFT+U **** 492 call psp_U_init() 493 494* **** initialize META GGA **** 495 call nwpw_meta_gga_init(control_gga()) 496 497* **** initialize vdw **** 498 call vdw_DF_init() 499 500 501* **** initialize pressure **** 502 calc_pressure = control_pressure().and.(control_version().eq.3) 503 pressure = 0.0d0 504 p1 = 0.0d0 505 p2 = 0.0d0 506 if (calc_pressure) then 507 call psp_stress_init() 508 call psp_stress_readall() 509 end if 510 511 512 513* |**************************| 514****************** summary of input data ********************** 515* |**************************| 516 call center_geom(cx,cy,cz) 517 call center_mass(gx,gy,gz) 518 call center_v_geom(vcx,vcy,vcz) 519 call center_v_mass(vgx,vgy,vgz) 520 mulliken = control_Mulliken() 521 522 if (taskid.eq.MASTER) then 523 write(luout,1111) np 524 write(luout,1117) np_i,np_j 525 if (mapping.eq.1) write(luout,1112) 526 if (mapping.eq.2) write(luout,1113) 527 if (mapping.eq.3) write(luout,1118) 528 if (control_balance()) then 529 write(luout,1114) 530 else 531 write(luout,1116) 532 end if 533 if (control_parallel_io()) then 534 write(luout,1119) 535 else 536 write(luout,1120) 537 end if 538 write(luout,1123) Parallel_maxthreads() 539 540 write(luout,1115) 541 write(luout,1121) control_boundry(),control_version() 542 if (ispin.eq.1) write(luout,1130) 'restricted' 543 if (ispin.eq.2) write(luout,1130) 'unrestricted' 544 545 call v_bwexc_print(luout,control_gga()) 546 547 if (fractional) write(luout,1132) 548 call pspw_print_SIC(luout) 549 call pspw_print_HFX(luout) 550 if (ion_makehmass2()) write(luout,1135) 551 write(luout,1140) 552 do ia = 1,ion_nkatm() 553 call psp_print(ia) 554c write(luout,1150) ia,ion_atom(ia), 555c > psp_zv(ia),psp_lmax(ia) 556c write(luout,1152) psp_lmax(ia) 557c write(luout,1153) psp_locp(ia) 558c write(luout,1154) psp_lmmax(ia) 559c if (control_version().eq.4) write(luout,1156) psp_rlocal(ia) 560c if (psp_semicore(ia)) 561c > write(luout,1155) psp_rcore(ia),psp_ncore(ia) 562c write(luout,1151) (psp_rc(i,ia),i=0,psp_lmax(ia)) 563 end do 564 565 icharge = -(ne(1)+ne(ispin)) 566 en(1) = ne(1) 567 en(ispin) = ne(ispin) 568 if (fractional) then 569 icharge = 0.0d0 570 do ms=1,ispin 571 en(ms) =0.0 572 do i=n1(ms),n2(ms) 573 icharge = icharge - (3-ispin)*dbl_mb(occ2(1)+i-1) 574 en(ms) = en(ms) + dbl_mb(occ2(1)+i-1) 575 end do 576 end do 577 end if 578 579 do ia=1,ion_nkatm() 580 icharge = icharge + ion_natm(ia)*psp_zv(ia) 581 end do 582 write(luout,1159) icharge 583 584 write(luout,1160) 585 write(luout,1170) (ion_atom(K),ion_natm(K),K=1,ion_nkatm()) 586 write(luout,1180) 587 do I=1,ion_nion() 588 if (ion_q_FixIon(I)) then 589 write(luout,1191) I,ion_aname(I), 590 > (ion_rion(K,I),K=1,3),ion_amass(i)/1822.89d0 591 else if (ion_q_xyzFixIon(I)) then 592 write(luout,1194) I,ion_aname(I),(ion_rion(K,I),K=1,3), 593 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 594 else 595 write(luout,1190) I,ion_aname(I), 596 > (ion_rion(K,I),K=1,3),ion_amass(i)/1822.89d0 597 end if 598 end do 599 write(luout,1200) cx,cy,cz 600 write(luout,1210) gx,gy,gz 601 602 603 write(luout,1181) 604 write(luout,1192) (I,ion_aname(I), 605 > (ion_vion(K,I),K=1,3),I=1,ion_nion()) 606 write(luout,1200) vcx,vcy,vcz 607 write(luout,1210) vgx,vgy,vgz 608 write(luout,1211) ion_nconstraints(),ion_ndof() 609 610 call pspw_charge_Print(luout) 611 call pspw_Efield_Print(luout) 612 613 if (fractional) then 614 write(luout,1219) en(1),en(ispin),' ( fractional)' 615 write(luout,1221) ne(1),neq(1), 616 > ne(ispin),neq(ispin),' (Fourier space)' 617 else 618 write(luout,1220) ne(1),neq(1), 619 > ne(ispin),neq(ispin),' (Fourier space)' 620 write(luout,1221) ne(1),neq(1), 621 > ne(ispin),neq(ispin),' (Fourier space)' 622 end if 623 write(luout,1230) 624 write(luout,1241) lattice_unita(1,1), 625 > lattice_unita(2,1), 626 > lattice_unita(3,1) 627 write(luout,1242) lattice_unita(1,2), 628 > lattice_unita(2,2), 629 > lattice_unita(3,2) 630 write(luout,1243) lattice_unita(1,3), 631 > lattice_unita(2,3), 632 > lattice_unita(3,3) 633 write(luout,1244) lattice_unitg(1,1), 634 > lattice_unitg(2,1), 635 > lattice_unitg(3,1) 636 write(luout,1245) lattice_unitg(1,2), 637 > lattice_unitg(2,2), 638 > lattice_unitg(3,2) 639 write(luout,1246) lattice_unitg(1,3), 640 > lattice_unitg(2,3), 641 > lattice_unitg(3,3) 642 write(luout,1231) lattice_omega() 643 write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 644 > pack_nwave_all(0),pack_nwave(0) 645 write(luout,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 646 > pack_nwave_all(1),pack_nwave(1) 647 if (control_version().eq.3) then 648 write(luout,1260) ewald_rcut(),ewald_ncut() 649 write(luout,1261) ewald_mandelung() 650 end if 651 652 if (control_has_ngrid_small()) then 653 write(luout,1229) 654 write(luout,1233) control_cell_name() 655 write(luout,1241) lattice_unita_small(1,1), 656 > lattice_unita_small(2,1), 657 > lattice_unita_small(3,1) 658 write(luout,1242) lattice_unita_small(1,2), 659 > lattice_unita_small(2,2), 660 > lattice_unita_small(3,2) 661 write(luout,1243) lattice_unita_small(1,3), 662 > lattice_unita_small(2,3), 663 > lattice_unita_small(3,3) 664 write(luout,1244) lattice_unitg_small(1,1), 665 > lattice_unitg_small(2,1), 666 > lattice_unitg_small(3,1) 667 write(luout,1245) lattice_unitg_small(1,2), 668 > lattice_unitg_small(2,2), 669 > lattice_unitg_small(3,2) 670 write(luout,1246) lattice_unitg_small(1,3), 671 > lattice_unitg_small(2,3), 672 > lattice_unitg_small(3,3) 673 call lattice_small_abc_abg(a,b,c,alpha,beta,gamma) 674 write(luout,1232) a,b,c,alpha,beta,gamma 675 write(luout,1231) lattice_omega_small() 676 write(luout,1250) lattice_ecut(), 677 > ngrid_small(1),ngrid_small(2),ngrid_small(3), 678 > pack_nwave_all(2),pack_nwave(2) 679 write(luout,1251) lattice_wcut(), 680 > ngrid_small(1),ngrid_small(2),ngrid_small(3), 681 > pack_nwave_all(3),pack_nwave(3) 682 end if 683 684 write(luout,1270) 685 if (.not.control_translation()) write(luout,1271) 686 if (.not.control_rotation()) write(luout,1272) 687 write(luout,1280) control_time_step(),control_fake_mass() 688 write(luout,1290) control_rte(),control_rti() 689 call ion_scaling_atoms_print(luout) 690 write(luout,1281) control_it_in()*control_it_out(), 691 > control_it_in(),control_it_out() 692 write(luout,1222) eke0,eki0,ekg 693 write(luout,1223) eke1,eki1 694 write(luout,1224) (eke1-eke0),(eki1-eki0) 695 if (control_Nose()) then 696 write(luout,1295) 697 if (Nose_restart()) then 698 write(luout,*) " thermostats resused" 699 else 700 write(luout,*) " thermostats initialized" 701 end if 702 do i=1,Nose_Mchain() 703 write(luout,1297) i,control_Nose_Te(),Nose_Qe(i), 704 > Nose_Pe(i),Nose_Ee0(i) 705 end do 706 do i=1,Nose_Nchain() 707 write(luout,1298) i,control_Nose_Tr(),Nose_Qr(i), 708 > Nose_Pr(i),Nose_Er0(i) 709 end do 710 else 711 write(luout,1294) 712 end if 713 if (calc_pressure) write(luout,1293) 714 if (control_SA()) then 715 write(luout,1296) sa_decay(1),sa_decay(2) 716 end if 717 718 719 if (mulliken) write(luout,1299) 720 write(luout,1300) 721 write(luout,1305) 722 call util_flush(luout) 723 end if 724 725* |***************************| 726****************** start iterations ********************** 727* |***************************| 728* **** open xyz and MOTION and dipole file **** 729 call xyz_init() ! unit=18 730 call MOTION_init(rtdb) ! unit=19 731 call dipole_motion_init(rtdb) ! unit=36 732 733* *** fei io **** 734 call fei_init(rtdb) 735 736 737* **** ecce print **** 738 call ecce_print_module_entry('task Car-Parrinello') 739 !call ecce_print_module_entry('driver') 740 call movecs_ecce_print_off() 741 742 743 744 745* ************************************ 746* **** open up other MOTION files **** 747* ************************************ 748 icount_shift = 0 749 750* **** open EMOTION file **** 751 if (.not.btdb_cget(rtdb,'cpmd:emotion_filename',1,filename)) 752 > call util_file_prefix('emotion',filename) 753 call util_file_name_noprefix(filename,.false., 754 > .false., 755 > full_filename) 756 if (taskid.eq.MASTER) then 757 758* **** check for backup file **** 759 call util_file_name_noprefix('EMOTION99-bak',.false., 760 > .false., 761 > full_bak) 762 inquire(file=full_bak,exist=found_bak) 763 if (found_bak) then 764 write(*,*) 765 write(*,*) "EMOTION99-bak exists:" 766 i=index(full_bak,' ') 767 j=index(full_filename,' ') 768 write(*,*) " Copying ",full_bak(1:i), 769 > " to ",full_filename(1:j) 770 write(*,*) 771 call util_file_copy(full_bak,full_filename) 772 end if 773 774 emotion_time_shift = 0.0d0 775 icount_shift = 0 776 inquire(file=full_filename,exist=found) 777 if (found) then 778 779* **** make a new backup file *** 780 call util_file_copy(full_filename,full_bak) 781 782 open(unit=31,file=full_filename,form='formatted', 783 > status='old') 784 do while (found) 785 read(31,*,end=100) emotion_time_shift,w,sum,gx,gy,gz 786 E(25) = E(25) + sum !*** take care of running sums *** 787 E(26) = E(26) + sum*sum 788 E(27) = E(27) + (sum+gx+gy) 789 E(28) = E(28) + (sum+gx+gy)**2 790 E(23) = E(23) + gz 791 E(24) = E(24) + gz*gz 792 icount_shift = icount_shift + 1 793 end do 794 100 continue 795#if defined(FUJITSU_SOLARIS) || defined(PSCALE) || defined(__crayx1) || defined(GCC46) 796 backspace 31 797#endif 798 else 799 open(unit=31,file=full_filename,form='formatted', 800 > status='new') 801 end if 802 end if 803 804 805* **** open EIGMOTION file **** 806 if (mulliken) then 807 if (.not.btdb_cget(rtdb,'cpmd:eigmotion_filename',1,filename)) 808 > call util_file_prefix('eigmotion',filename) 809 call util_file_name_noprefix(filename,.false., 810 > .false., 811 > full_filename) 812 if (taskid.eq.MASTER) 813 > open(unit=32,file=full_filename,form='formatted') 814 end if 815 816* **** open HMOTION file **** 817 if (mulliken) then 818 if (.not.btdb_cget(rtdb,'cpmd:hmotion_filename',1,filename)) 819 > call util_file_prefix('hmotion',filename) 820 call util_file_name_noprefix(filename,.false., 821 > .false., 822 > full_filename) 823 if (taskid.eq.MASTER) 824 > open(unit=34,file=full_filename,form='formatted') 825 end if 826 827* **** open OMOTION file **** 828 if (mulliken) call Orb_Init(rtdb,ispin,ne) !unit=33 829 830c* **** write initial position to xyz data **** 831c call xyz_write() 832 833* ***** first step using velocity **** 834 verlet = .false. 835 call inner_loop_md(verlet,sa_alpha,ispin,ne,neq, 836 > npack1,nfft3d,nemaxq, 837 > dcpl_mb(psi0(1)), 838 > dcpl_mb(psi1(1)), 839 > dcpl_mb(psi2(1)), 840 > dbl_mb(dn(1)), 841 > 1,0,E, 842 > dbl_mb(hml(1)),dbl_mb(lmd(1)),dbl_mb(lmd1(1)), 843 > dcpl_mb(Hpsi(1)),dcpl_mb(psir(1)), 844 > calc_pressure,pressure,p1,p2, 845 > fractional, 846 > dbl_mb(occ0(1)),dbl_mb(occ1(1)),dbl_mb(occ2(1))) 847 848 849 if (taskid.eq.MASTER) call current_second(cpu2) 850 if ((taskid.eq.MASTER).and.(.not.calc_pressure)) 851 > CALL nwpw_message(6) 852 if ((taskid.eq.MASTER).and.(calc_pressure)) 853 > CALL nwpw_message(9) 854 855 it_in = control_it_in() 856 it_out = control_it_out() 857 icount = 0 858 verlet = .true. 859 eke = 0.0d0 860 dt = control_time_step() 861 862 if (it_out.lt.1) goto 102 863 864 Te_new = Te_init 865 Tr_new = Tr_init 866 101 continue 867 icount = icount + 1 868 call inner_loop_md(verlet,sa_alpha,ispin,ne,neq, 869 > npack1,nfft3d,nemaxq, 870 > dcpl_mb(psi0(1)), 871 > dcpl_mb(psi1(1)), 872 > dcpl_mb(psi2(1)), 873 > dbl_mb(dn(1)), 874 > it_in,((icount-1)*it_in), 875 > E, 876 > dbl_mb(hml(1)),dbl_mb(lmd(1)),dbl_mb(lmd1(1)), 877 > dcpl_mb(Hpsi(1)), dcpl_mb(psir(1)), 878 > calc_pressure,pressure,p1,p2, 879 > fractional, 880 > dbl_mb(occ0(1)),dbl_mb(occ1(1)),dbl_mb(occ2(1))) 881 882 eke = eke + E(3) 883 884 885 !**** metadynamics and tamd update **** 886 call meta_update(ispin,neq,dcpl_mb(psi1(1)),E) 887 call tamd_update(ispin,neq,dcpl_mb(psi1(1)),E) 888 889 if (taskid.eq.MASTER) then 890 891 if (calc_pressure) then 892 if (SA) then 893 write(luout,1309) icount*it_in,E(1),E(2),E(3),E(4), 894 > Te_new,Tr_new,pressure*autoatm 895 else 896 write(luout,1310) icount*it_in,E(1),E(2),E(3),E(4), 897 > ion_Temperature(),pressure*autoatm 898 end if 899 else 900 if (SA) then 901 write(luout,1309) icount*it_in,E(1),E(2),E(3),E(4), 902 > Te_new,Tr_new 903 else 904 write(luout,1310) icount*it_in,E(1),E(2),E(3),E(4), 905 > ion_Temperature() 906 end if 907 end if 908 call util_flush(luout) 909 910* **** write out EMOTION data **** 911 qave = E(23)/dble(icount+icount_shift) 912 qvar = E(24)/dble(icount+icount_shift) 913 qvar = qvar - qave*qave 914 eave = E(25)/dble(icount+icount_shift) 915 evar = E(26)/dble(icount+icount_shift) 916 evar = evar - eave*eave 917 have = E(27)/dble(icount+icount_shift) 918 hvar = E(28)/dble(icount+icount_shift) 919 hvar = hvar - have*have 920 if (control_Nose()) then 921 write(31,1311) icount*it_in*dt + emotion_time_shift, 922 > e(1),e(2),e(3),e(4),e(14),e(5),e(6), 923 > e(7),e(8),e(9),e(10), 924 > eave,evar,have,hvar,qave,qvar, 925 > ion_Temperature(),pressure 926 else 927 write(31,1311) icount*it_in*dt + emotion_time_shift, 928 > e(1),e(2),e(3),e(4),e(14),e(5),e(6), 929 > e(7),e(8), 930 > eave,evar,have,hvar,qave,qvar, 931 > ion_Temperature(),pressure 932 end if 933 call util_flush(31) 934 935* **** write out EIGMOTION data -diagonal hml matrix **** 936 if (mulliken) then 937 write(32,1311) icount*it_in*dt, 938 > (( dbl_mb(hml(1)+ii-1+(ii-1)*ne(1)+(ms-1)*ne(1)*ne(1)), 939 > ii=1,ne(ms)),ms=1,ispin) 940 call util_flush(32) 941 end if 942 943* **** write out HMOTION data - hml matrix **** 944 if (mulliken) then 945 write(34,1312) icount*it_in*dt,ispin 946 do ms=1,ispin 947 write(34,1313) ms,ne(ms),ne(ms) 948 do ii=1,ne(ms) 949 write(34,1311) 950 > (dbl_mb(hml(1)+ii-1+(jj-1)*ne(1)+(ms-1)*ne(1)*ne(1)), 951 > jj=1,ne(ms)) 952 end do 953 end do 954 call util_flush(34) 955 end if 956 957 end if 958 959 960* **** write xyz, MOTION, and dipole data **** 961 call xyz_write() 962 call MOTION_write(icount*it_in*control_time_step()) 963 964 call dipole_motion_write((control_version().eq.3), 965 > (icount*it_in*control_time_step()), 966 > ispin,ne,neq,npack1,nfft3d, 967 > dbl_mb(dn(1)), 968 > dcpl_mb(psi1(1))) 969 970 971* **** write OMOTION data **** 972 if (mulliken) call Orb_Write(dcpl_mb(psi1(1))) 973 974* **** update thermostats using SA decay **** 975 if (SA) then 976 if(abs(sa_decay(1)).gt.tollz) 977 * t1 = icount*it_in*dt/sa_decay(1) 978 if(abs(sa_decay(2)).gt.tollz) 979 * t2 = icount*it_in*dt/sa_decay(2) 980 Te_new = Te_init*dexp(-t1) 981 Tr_new = Tr_init*dexp(-t2) 982 call Nose_reset_T(Te_new,Tr_new) 983 end if 984 985 986* **** exit early **** 987 if (control_out_of_time()) then 988 if (taskid.eq.MASTER) 989 > write(luout,*) ' *** out of time. iteration terminated' 990 go to 102 991 end if 992 if (icount.lt.it_out) go to 101 993 if (taskid.eq.MASTER) 994 > write(luout,*) 995 > '*** arrived at the Maximum iteration. terminated.' 996 997*:::::::::::::::::::: end of iteration loop ::::::::::::::::::::::::: 998 999 102 continue 1000 1001* **** close xyz,MOTION and dipole files **** 1002 call xyz_end() 1003 call MOTION_end() 1004 call dipole_motion_end() 1005 if (taskid.eq.MASTER) then 1006 close(unit=31) 1007 close(unit=32) 1008 close(unit=34) 1009 1010* **** remove EMOTION backup file *** 1011 call util_file_name_noprefix('EMOTION99-bak',.false., 1012 > .false., 1013 > full_bak) 1014 call util_file_unlink(full_bak) 1015 end if 1016 1017* *** close fei io **** 1018 call fei_end() 1019 1020* **** close OMOTION file **** 1021 if (mulliken) call Orb_End() 1022 1023* **** ecce print **** 1024 !call ecce_print_module_exit('driver', 'ok') 1025 call ecce_print_module_exit('task Car-Parrinello', 'ok') 1026 1027 1028* **** finalize pressure **** 1029 if (calc_pressure) then 1030 call psp_stress_end() 1031 end if 1032 1033 1034 if (taskid.eq.MASTER) CALL nwpw_message(3) 1035 if (taskid.eq.MASTER) call current_second(cpu3) 1036 1037 1038* |****************************************| 1039*********** produce CHECK file and diagonalize hml ***************** 1040* |****************************************| 1041 1042* **** produce CHECK FILE **** 1043 if (taskid.eq.MASTER) then 1044 call util_file_name('CHECK',.true., 1045 > .false., 1046 > full_filename) 1047 open(unit=17,file=full_filename,form='formatted') 1048 end if 1049 1050* **** check total number of electrons **** 1051 do ms =1,ispin 1052 call D3dB_r_dsum(1,dbl_mb(dn(1)+(ms-1)*n2ft3d),sumall) 1053 en1(ms) = sumall*lattice_omega() 1054 > /dble(ngrid(1)*ngrid(2)*ngrid(3)) 1055 end do 1056 1057 if (psp_pawexist()) then 1058 if (.not.BA_push_get(mt_dbl,n2ft3d,'tmp1',tmp1(2),tmp1(1))) 1059 > call errquit( 1060 > 'cpmdv5: out of stack memory',0,MA_ERR) 1061 1062 call psp_qlm_atom(ispin,neq,dcpl_mb(psi1(1))) 1063 do ms=1,ispin 1064 call nwpw_compcharge_gen_dn_cmp_smooth_ms(ms,dbl_mb(tmp1(1))) 1065 call Pack_c_unpack(0,dbl_mb(tmp1(1))) 1066 call D3dB_cr_fft3b(1,dbl_mb(tmp1(1))) 1067 call D3dB_r_Zero_Ends(1,dbl_mb(tmp1(1))) 1068 call D3dB_r_dsum(1,dbl_mb(tmp1(1)),sumall) 1069 en2(ms) = sumall*lattice_omega() 1070 > /dble(ngrid(1)*ngrid(2)*ngrid(3)) 1071 end do 1072 if (.not.BA_pop_stack(tmp1(2))) 1073 > call errquit( 1074 > 'cpmdv5: popping stack memory',0,MA_ERR) 1075 else 1076 en2(1) = 0.0d0 1077 en2(2) = 0.0d0 1078 end if 1079 en(1) = en1(1)+en2(1) 1080 en(2) = en1(2)+en2(2) 1081 1082 if (taskid.eq.MASTER) then 1083 write(17,1320) (en(ms),ms=1,ispin) 1084 if (psp_pawexist()) then 1085 write(17,1322) (en1(ms),ms=1,ispin) 1086 write(17,1323) (en2(ms),ms=1,ispin) 1087 end if 1088 end if 1089 1090* **** comparison between hamiltonian an lambda matrix **** 1091 if (taskid.eq.MASTER) write(17,1330) 1092 do ms=1,ispin 1093 do i=1,ne(ms) 1094 do j=1,ne(ms) 1095 w = Dneall_m_value(0,ms,i,j,dbl_mb(hml(1))) 1096 sum = Dneall_m_value(0,ms,i,j,dbl_mb(lmd(1))) 1097 1098 if (taskid.eq.MASTER) 1099 > write(17,1340) ms,i,j,w,sum,w-sum 1100 1101 end do 1102 end do 1103 end do 1104 1105* **** check orthonormality **** 1106 if (taskid.eq.MASTER) then 1107 write(17,1350) 1108 end if 1109 1110 call Dneall_ffm_Multiply(0,dcpl_mb(psi1(1)), 1111 > dcpl_mb(psi1(1)),npack1, 1112 > dbl_mb(lmd(1))) 1113 do ms=1,ispin 1114 do j=1,ne(ms) 1115 do i=j,ne(ms) 1116 w = Dneall_m_value(0,ms,i,j,dbl_mb(lmd(1))) 1117 if (taskid.eq.MASTER) write(17,1360) ms,i,j,w 1118 end do 1119 end do 1120 end do 1121 1122* **** close check file **** 1123 if (taskid.eq.MASTER) then 1124 close(17) 1125 end if 1126 1127 1128 1129* ***** do not diagonalize the hamiltonian matrix ***** 1130 if (pspw_SIC()) then 1131 call dcopy(2*npack1*nemaxq, 1132 > dcpl_mb(psi1(1)),1, 1133 > dcpl_mb(psi2(1)),1) 1134 1135* ***** diagonalize the hamiltonian matrix ***** 1136 else 1137 1138c if (fractional) then 1139c call Dneall_m_HmltimesSA(0,dbl_mb(hml(1)),dbl_mb(fweight(1))) 1140c end if 1141 1142 call Dneall_m_diagonalize(0,dbl_mb(hml(1)), 1143 > dbl_mb(eig(1)),.false.) 1144 1145 1146c if (fractional) then 1147c do ii=1,ne(ms) 1148c dbl_mb(eig(1)+(ii-1)+(ms-1)*n) 1149c > =dbl_mb(eig(1)+(ii-1)+(ms-1)*n) 1150c > /dbl_mb(fweight(1)+(ii-1)+(ms-1)*n) 1151c end do 1152c end if 1153 1154 1155* **** do not rotate for wannier localization algorithm **** 1156 if (.not.pspw_HFX_localize2()) then 1157* *** rotate current psi *** 1158 call Dneall_fmf_Multiply(0,dcpl_mb(psi1(1)),npack1, 1159 > dbl_mb(hml(1)), 1.0d0, 1160 > dcpl_mb(psi2(1)),0.0d0) 1161 1162 1163* *** rotate current v_psi *** 1164 call dcopy(2*npack1*nemaxq,dcpl_mb(psi0(1)),1, 1165 > dcpl_mb(psi1(1)),1) 1166 1167 call Dneall_fmf_Multiply(0,dcpl_mb(psi1(1)),npack1, 1168 > dbl_mb(hml(1)), 1.0d0, 1169 > dcpl_mb(psi0(1)),0.0d0) 1170 end if 1171 1172 end if 1173 1174 1175 1176* |***************************| 1177****************** report summary of results ********************** 1178* |***************************| 1179 call center_geom(cx,cy,cz) 1180 call center_mass(gx,gy,gz) 1181 call center_v_geom(vcx,vcy,vcz) 1182 call center_v_mass(vgx,vgy,vgz) 1183 1184 if (taskid.eq.MASTER) then 1185 call print_elapsed_time(icount*it_in*dt) 1186 write(luout,1300) 1187 write(luout,1410) 1188 write(luout,1420) 1189 do I=1,ion_nion() 1190 if (ion_q_FixIon(I)) then 1191 write(luout,1191) I,ion_aname(I),(ion_rion(k,i),K=1,3), 1192 > ion_amass(I)/1822.89d0 1193 else if (ion_q_xyzFixIon(I)) then 1194 write(6,1194) I,ion_aname(I),(ion_rion(k,i),K=1,3), 1195 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 1196 else 1197 write(luout,1190) I,ion_aname(I),(ion_rion(k,i),K=1,3), 1198 > ion_amass(I)/1822.89d0 1199 end if 1200 end do 1201 write(luout,1200) cx,cy,cz 1202 write(luout,1210) gx,gy,gz 1203 1204 1205 write(luout,1421) 1206 write(luout,1192) (I,ion_aname(I), 1207 > (ion_vion(K,I),K=1,3),I=1,ion_nion()) 1208 write(luout,1200) vcx,vcy,vcz 1209 write(luout,1210) vgx,vgy,vgz 1210 write(luout,1211) ion_nconstraints(),ion_ndof() 1211 1212 call pspw_charge_Print(luout) 1213 call pspw_Efield_Print(luout) 1214 1215 write(luout,*) 1216 write(luout,1320) en(1),en(ispin),' (real space)' 1217 if (psp_pawexist()) then 1218 write(luout,1322) en1(1),en1(ispin),' (real space)' 1219 write(luout,1323) en2(1),en2(ispin),' (real space)' 1220 end if 1221 1222 if (psp_pawexist()) then 1223 write(luout,1434) (E(1)+E(36)+E(45)), 1224 > (E(1)+E(36)+E(45))/ion_nion() 1225 end if 1226 1227* **** write APC potential and charges *** 1228 if (pspw_V_APC_on()) call pspw_shortprint_APC(luout) 1229 1230 write(luout,1430) E(2),E(2)/ion_nion() 1231 write(luout,1440) E(5),E(5)/n2(ispin) 1232 write(luout,1450) E(6),E(6)/n2(ispin) 1233 write(luout,1460) E(7),E(7)/n2(ispin) 1234 if (pspw_SIC()) then 1235 write(luout,1455) E(16),E(16)/n2(ispin) 1236 write(luout,1456) E(17),E(17)/n2(ispin) 1237 end if 1238 if (pspw_HFX()) then 1239 write(luout,1457) E(20),E(20)/n2(ispin) 1240 end if 1241 if (psp_U_psputerm()) then 1242 write(luout,1458) E(29),E(29)/n2(ispin) 1243 end if 1244 if (meta_found()) then 1245 write(luout,1459) E(31),E(31)/ion_nion() 1246 end if 1247 if (tamd_found()) then 1248 write(luout,1461) E(34),E(34)/ion_nion() 1249 end if 1250 if (pspw_V_APC_on()) then 1251 write(luout,1505) E(52),E(52)/ion_nion() 1252 end if 1253 write(luout,1470) E(8),E(8)/ion_nion() 1254 write(luout,1471) E(3),E(3)/n2(ispin) 1255 write(luout,1472) ion_ke(),ion_ke()/ion_nion() 1256 1257 1258 if (pspw_qmmm_found()) then 1259 write(luout,1700) 1260 write(luout,1701) 1261 write(luout,1702) E(11) 1262 write(luout,1703) E(12) 1263 write(luout,1704) E(13) 1264 qave = E(23)/dble(icount+icount_shift) 1265 qvar = E(24)/dble(icount+icount_shift) 1266 qvar = qvar - qave*qave 1267 write(luout,1707) pspw_qmmm_lambda() 1268 write(luout,1705) E(14),qave,qvar 1269 !write(luout,1706) qave,qvar 1270 end if 1271 if (ion_disp_on()) then 1272 write(luout,1720) E(33) 1273 end if 1274 1275 if (field_exist) then 1276 write(luout,1800) 1277 write(luout,1801) 1278 write(luout,1805) E(19)+E(20)+E(21) 1279 write(luout,1802) E(19) 1280 write(luout,1803) E(20) 1281 write(luout,1804) E(21) 1282 end if 1283 1284 if (control_Nose()) then 1285 write(luout,1473) E(9),E(9)/n2(ispin) 1286 write(luout,1474) E(10),E(10)/ion_nion() 1287 end if 1288 write(luout,1226) E(3),ion_ke(),ion_com_ke() 1289 eke = eke/dble(it_out) 1290 eke = 2.0d0*eke/kb/(ne(1)+ne(ispin))/pack_nwave_all(1) 1291 !eke = 2.0d0*eke/kb/(ne(1)+ne(ispin)) 1292 1293* **** write out Temperatures **** 1294 write(luout,1491) eke 1295 eki0 = ion_Temperature() 1296c if (pspw_qmmm_found()) then 1297c eki1 =pspw_qmmm_Temperature() 1298c sum = ion_nion() + pspw_qmmm_nion() - 2.0d0 1299c eki0 = eki0*((ion_nion()-2.0d0)/sum) 1300c > + eki1*((pspw_qmmm_nion()-2.0d0)/sum) 1301c end if 1302 write(luout,1480) eki0 1303 write(luout,1490) ion_com_Temperature() 1304 1305 eave = E(25)/dble(icount+icount_shift) 1306 evar = E(26)/dble(icount+icount_shift) 1307 evar = evar - eave*eave 1308 have = E(27)/dble(icount+icount_shift) 1309 hvar = E(28)/dble(icount+icount_shift) 1310 hvar = hvar - have*have 1311 cv = (evar)/(kb*ion_Temperature()**2) 1312 cv = cv/dble(ion_nion()) 1313 write(luout,1492) eave,have 1314 write(luout,1493) evar,hvar 1315 write(luout,1494) cv 1316 1317* **** write out diagonal <psi|H|psi> matrix **** 1318 if (pspw_SIC()) then 1319 1320 n = ne(1) 1321 nn = n*n 1322 do ms=1,ispin 1323 if (ms.eq.1) write(luout,1331) 1324 if (ms.eq.2) write(luout,1332) 1325 !*** call Gainsville matrix output *** 1326 call output(dbl_mb(hml(1)+(ms-1)*nn), 1327 > 1,ne(ms),1,ne(ms), 1328 > n,n,1) 1329 end do 1330 1331 1332* **** write out KS eigenvalues **** 1333 else 1334 write(luout,1500) 1335 NN=NE(1)-NE(2) 1336 EV=27.2116d0 1337 1338 if (fractional) then 1339 do i=1,NN 1340 write(luout,1511) dbl_mb(EIG(1)+i-1), 1341 > dbl_mb(EIG(1)+i-1)*EV, 1342 > dbl_mb(occ2(1)+i-1) 1343 end do 1344 do i=1,ne(2) 1345 write(luout,1511) dbl_mb(EIG(1)+i-1+NN), 1346 > dbl_mb(EIG(1)+i-1+NN)*EV, 1347 > dbl_mb(occ2(1)+i-1+NN), 1348 > dbl_mb(EIG(1)+i-1+n1(2)-1), 1349 > dbl_mb(EIG(1)+i-1+n1(2)-1)*EV, 1350 > dbl_mb(occ2(1)+i-1+n1(2)-1) 1351 end do 1352 else 1353 do i=1,NN 1354 write(luout,1510) dbl_mb(EIG(1)+i-1),dbl_mb(EIG(1)+i-1)*EV 1355 end do 1356 do i=1,ne(2) 1357 write(luout,1510) dbl_mb(EIG(1)+i-1+NN), 1358 > dbl_mb(EIG(1)+i-1+NN)*EV, 1359 > dbl_mb(EIG(1)+i-1+n1(2)-1), 1360 > dbl_mb(EIG(1)+i-1+n1(2)-1)*EV 1361 end do 1362 end if 1363 1364 end if 1365 end if 1366 1367* **** write out extended Born solvation energies **** 1368 if (nwpw_born_on()) then 1369 Egas = control_gas_energy() 1370 if (taskid.eq.MASTER) then 1371 write(luout,1740) 1372 write(luout,1741) nwpw_born_screen() 1373 write(luout,1745) E(52),E(52)*27.2116d0*23.06d0 1374 if (dabs(Egas).gt.1.0d-6) 1375 > write(luout,1746) E(2)-Egas, 1376 > (E(2)-Egas)*27.2116d0*23.06d0 1377 call nwpw_born_print(luout,Egas,E(2)) 1378 end if 1379 end if 1380 1381 if (taskid.eq.MASTER) then 1382* *** Extra energy output added for QA test **** 1383 write(luout,1600) E(2) 1384 end if 1385 1386* |***************************| 1387****************** Prologue ********************** 1388* |***************************| 1389 1390* **** calculate spin contamination **** 1391 call Calculate_psi_spin2(ispin,ne,npack1,dcpl_mb(psi2(1)), 1392 > fractional,dbl_mb(occ2(1)), 1393 > dipole) 1394 1395* **** calculate the Dipole *** 1396 call Calculate_Dipole(ispin,ne,n2ft3d,dbl_mb(dn(1)),dipole) 1397 1398 1399* ***** write wavefunctions and v_wavefunctions **** 1400 call psi_write(ispin,ne,dcpl_mb(psi2(1)), 1401 > smearoccupation,dbl_mb(occ2(1))) 1402 call v_psi_write(ispin,ne,dcpl_mb(psi0(1))) 1403 1404* **** write geometry to rtdb **** 1405 call pspw_charge_write(rtdb) 1406 call ion_write(rtdb) 1407 1408* **** deallocate heap memory **** 1409 if (control_version().eq.3) call ewald_end() 1410 call strfac_end() 1411 if (controL_version().eq.3) call coulomb_end() 1412 if (controL_version().eq.4) call coulomb2_end() 1413 call ke_end() 1414 call mask_end() 1415 call Pack_end() 1416 call G_end() 1417 call psp_U_end() 1418 call vdw_DF_end() 1419 call nwpw_meta_gga_end() 1420 call pspw_end_SIC() 1421 call pspw_end_HFX() 1422 call pspw_end_APC() 1423 call pspw_qmmm_end() 1424 call meta_finalize(rtdb) 1425 call tamd_finalize(rtdb) 1426 call dplot_iteration_end() 1427 call pspw_charge_end() 1428 call pspw_Efield_end() 1429c call frac_occ_end() 1430 if (control_Nose()) call Nose_end() 1431 if (psp_pawexist()) call nwpw_gintegrals_end() 1432 !if (psp_pawexist()) call nwpw_cgintegrals_end() 1433 1434 call ion_end() 1435 call psp_end() 1436 call ion_end_FixIon() 1437 1438 value = BA_free_heap(psir(2)) 1439 value = BA_free_heap(Hpsi(2)) 1440 value = BA_free_heap(dn(2)) 1441 value = BA_free_heap(eig(2)) 1442 value = Dneall_m_free(hml) 1443 value = Dneall_m_free(lmd) 1444 value = Dneall_m_free(lmd1) 1445 value = BA_free_heap(psi0(2)) 1446 value = BA_free_heap(psi1(2)) 1447 value = BA_free_heap(psi2(2)) 1448 if (fractional) then 1449 value = BA_free_heap(occ0(2)) 1450 value = BA_free_heap(occ1(2)) 1451 value = BA_free_heap(occ2(2)) 1452 end if 1453 call D3dB_pfft_end() 1454 call D3dB_end(1) 1455 if (control_version().eq.4) call D3dB_end(2) 1456 if (control_has_ngrid_small()) call D3dB_end(3) 1457 call Dne_end() 1458 call psi_data_end() 1459 1460* **** do anaylysis on MOTION files **** 1461 call cpmd_properties(rtdb) 1462 1463 1464* |***************************| 1465****************** report consumed cputime ********************** 1466* |***************************| 1467 if (taskid.eq.MASTER) then 1468 CALL current_second(cpu4) 1469 1470 T1=CPU2-CPU1 1471 T2=CPU3-CPU2 1472 T3=CPU4-CPU3 1473 T4=CPU4-CPU1 1474 AV=T2/dble(icount*it_in) 1475 write(luout,*) 1476 write(luout,*) '-----------------' 1477 write(luout,*) 'cputime in seconds' 1478 write(luout,*) 'prologue : ',T1 1479 write(luout,*) 'main loop : ',T2 1480 write(luout,*) 'epilogue : ',T3 1481 write(luout,*) 'total : ',T4 1482 write(luout,*) 'cputime/step: ',AV 1483 write(luout,*) 1484 1485 call nwpw_timing_print_final(.true.,(icount*it_in)) 1486 CALL nwpw_message(4) 1487 end if 1488 1489 1490 call Parallel2d_Finalize() 1491 call Parallel_Finalize() 1492 cpmdv5 = value 1493 return 1494 1495 1496*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 1497 1000 FORMAT(10X,'****************************************************') 1498 1010 FORMAT(10X,'* *') 1499 1020 FORMAT(10X,'* Car-Parrinello microcluster calculation *') 1500 1030 FORMAT(10X,'* [ extended Lagrangian molecular ] *') 1501 1031 FORMAT(10X,'* [ dynamics simulation ] *') 1502 1035 FORMAT(10x,'* [ NorthWest Chemistry implementation ] *') 1503 1040 FORMAT(10X,'* version #5.00 06/01/99 *') 1504 1041 FORMAT(10X,'* This code was developed by Eric J. Bylaska, *') 1505 1042 FORMAT(10X,'* and was based upon algorithms and code *') 1506 1043 FORMAT(10X,'* developed by the group of Prof. John H. Weare *') 1507 1100 FORMAT(//) 1508 1110 FORMAT(10X,'================ input data ========================') 1509 1111 FORMAT(/' number of processors used:',I10) 1510 1112 FORMAT( ' parallel mapping : 1d-slab') 1511 1113 FORMAT( ' parallel mapping : 2d-hilbert') 1512 1114 FORMAT( ' parallel mapping : balanced') 1513 1115 FORMAT(/' options:') 1514 1116 FORMAT( ' parallel mapping : not balanced') 1515 1117 FORMAT( ' processor grid :',I4,' x',I4) 1516 1118 FORMAT( ' parallel mapping : 2d-hcurve') 1517 1119 FORMAT( ' parallel io : on') 1518 1120 FORMAT( ' parallel io : off') 1519 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 1520 1123 FORMAT( ' number of threads :',I10) 1521 1130 FORMAT(5X,' electron spin = ',A) 1522 1131 FORMAT(5X,' exchange-correlation = ',A) 1523 1132 FORMAT(5X,' using fractional occupation') 1524c 1135 FORMAT(/' The masses of QM H atoms converted to 2.0 amu. ', 1525c > /' To turn off this default', 1526c > ' set nwpw:makehmass2 .false.') 1527 1135 FORMAT(/' The masses of QM H atoms converted to 2.0 amu. ', 1528 > /' To turn off this default', 1529 > /' nwpw', 1530 > /' makehmass2 off', 1531 > /' end') 1532 1140 FORMAT(/' elements involved in the cluster:') 1533 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I1) 1534 1151 FORMAT(5X,' cutoff =',4F8.3) 1535 1152 FORMAT(12X,' highest angular component : ',i2) 1536 1153 FORMAT(12X,' local potential used : ',i2) 1537 1154 FORMAT(12X,' number of non-local projections: ',i2) 1538 1155 FORMAT(12X,' semicore corrections included : ', 1539 > F6.3,' (radius) ',F6.3,' (charge)') 1540 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 1541 1159 FORMAT(/' total charge=',F8.3) 1542 1160 FORMAT(/' atomic composition:') 1543 1170 FORMAT(7(5X,A4,':',I5)) 1544 1180 FORMAT(/' initial position of ions:') 1545 1181 FORMAT(/' initial velocity of ions:') 1546 1190 FORMAT(5X, I4, A5 ,' (',3F11.5,' ) - atomic mass= ',F7.3,' ') 1547 1191 FORMAT(5X, I4, A5 ,' (',3F11.5, 1548 > ' ) - atomic mass= ',F7.3,' - fixed') 1549 1192 FORMAT(5X, I4, A5 ,' (',3F11.5,' )') 1550 1193 FORMAT(5X, I4, A5, ' (',3F11.5, 1551 > ' ) - atomic mass= ',F7.3,' - z fixed') 1552 1194 FORMAT(5X, I4, A5, ' (',3F11.5, 1553 > ' ) - atomic mass= ',F7.3,A) 1554 1555 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 1556 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 1557 1211 FORMAT(5X,' number of constraints = ', I6,' ( DOF = ',I6,' )' ) 1558 1219 FORMAT(/' number of electrons: spin up=',F6.2,' down=',F6.2,A) 1559 1220 FORMAT(/' number of electrons: spin up=',I6, 1560 > ' (',I4,' per task)', 1561 > ' down=',I6, 1562 > ' (',I4,' per task)', 1563 > A) 1564 1221 FORMAT( ' number of orbitals : spin up=',I6, 1565 > ' (',I4,' per task)', 1566 > ' down=',I6, 1567 > ' (',I4,' per task)', 1568 > A) 1569 1222 format(5x,' initial kinetic energy: ',e12.5,' (psi)', 2x, 1570 > e12.5,' (ion)',/50x, 1571 > e12.5,' (c.o.m.)') 1572 1223 format(5x,' after scaling: ',e12.5,' (psi)', 2x, 1573 > e12.5,' (ion)') 1574 1224 format(5x,' increased energy: ',e12.5,' (psi)', 2x, 1575 > e12.5,' (ion)') 1576 1226 format(/' final kinetic energy: ',e12.5,' (psi)', 2x, 1577 > e12.5,' (ion)',/44x, 1578 > e12.5,' (c.o.m.)') 1579 1229 FORMAT(/' small supercell:') 1580 1230 FORMAT(/' supercell:') 1581 1231 FORMAT(5x,' volume : ',F12.1) 1582 1232 FORMAT(5x,' lattice: a= ',f8.3,' b= ',f8.3,' c= ',f8.3, 1583 > /5x,' alpha=',f8.3,' beta=',f8.3,' gamma=',f8.3) 1584 1233 FORMAT(5x,' cell_name: ',A) 1585 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 1586 1242 FORMAT(5x,' a2=<',3f8.3,' >') 1587 1243 FORMAT(5x,' a3=<',3f8.3,' >') 1588 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >') 1589 1245 FORMAT(5x,' b2=<',3f8.3,' >') 1590 1246 FORMAT(5x,' b3=<',3f8.3,' >') 1591 1592 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1593 & '( ',I8,' waves ',I8,' per task)') 1594 1251 FORMAT(5X,' wavefnc cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1595 & '( ',I8,' waves ',I8,' per task)') 1596 1597 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,' and',I3) 1598 1261 FORMAT(5X,' madelung=',f14.8) 1599 1270 FORMAT(/' technical parameters:') 1600 1271 FORMAT(5x, ' translation constrained') 1601 1272 FORMAT(5x, ' rotation constrained') 1602 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1) 1603 1281 FORMAT(5X, ' maximum iterations =',I10, 1604 > ' ( ',I4,' inner ',I6,' outer )') 1605 1290 FORMAT(5X, ' cooling/heatting rates: ',e12.5,' (psi)',2x, 1606 > e12.5,' (ion)') 1607 1608 1293 format(/' Pressure Output Generated ') 1609 1294 format(/' Constant Energy Simulation ') 1610 1295 format(/' Nose-Hoover Simulation - Thermostat Parameters:') 1611 1296 format(5x, 'SA decay rates =',e10.3,' (elc)',e10.3,' (ion)') 1612 1297 format(5x, 'link = ',I3, 1613 > ' Te =',f8.2,' Qe =',e10.3,' 2*pi/we=',e10.3,' Ee0=',e10.3) 1614 1298 format(5x, 'link = ',I3, 1615 > ' Tr =',f8.2,' Qr =',e10.3,' 2*pi/wr=',e10.3,' Er0=',e10.3) 1616 1299 format(//' Mulliken Analysis Output Generated ') 1617 1300 FORMAT(//) 1618 1305 FORMAT(10X,'============ Car-Parrinello iteration ==============') 1619 1309 FORMAT(I8,2E19.10,2E14.5,2F9.1,3E11.3) 1620 1310 FORMAT(I8,2E19.10,2E14.5,F14.2,3E11.3) 1621 1311 format(100e19.10) 1622 1312 format(e14.6,i3) 1623 1313 format(3i4) 1624 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 1625 1321 FORMAT(' total charge of system:',F11.5,A) 1626 1322 FORMAT(' plane-wave part: ',F11.5,' ',F11.5,A) 1627 1323 FORMAT(' augmented part: ',F11.5,' ',F11.5,A) 1628 1330 FORMAT(/' comparison between hamiltonian and lambda matrix') 1629 1331 FORMAT(/' Elements of Hamiltonian matrix (up/restricted)') 1630 1332 FORMAT(/' Elements of Hamiltonian matrix (down)') 1631 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7) 1632 1341 FORMAT(I3,2I3,' H=',E16.6) 1633 1350 FORMAT(/' orthonormality') 1634 1360 FORMAT(I3,2I3,E18.7) 1635 1370 FORMAT(I3) 1636 1380 FORMAT(' ''',a,'''',I4) 1637 1390 FORMAT(I3) 1638 1400 FORMAT(I3,3E18.8/3X,3E18.8) 1639 1410 FORMAT(10X,'============= summary of results =================') 1640 1420 FORMAT(/' final position of ions:') 1641 1421 FORMAT(/' final velocity of ions:') 1642 1430 FORMAT(/' total energy :',E19.10,' (',E15.5,'/ion)') 1643 1431 FORMAT(/' QM Energies') 1644 1432 FORMAT( '------------') 1645 1434 FORMAT(//' total paw energy :',E19.10,' (',E15.5,'/ion)') 1646 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)') 1647 1450 FORMAT( ' hartree energy :',E19.10,' (',E15.5,'/electron)') 1648 1455 FORMAT( ' SIC-hartree energy :',E19.10,' (',E15.5,'/electron)') 1649 1456 FORMAT( ' SIC-exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1650 1457 FORMAT( ' HF exchange energy :',E19.10,' (',E15.5,'/electron)') 1651 1458 FORMAT( ' DFT+U energy :',E19.10,' (',E15.5,'/ion)') 1652 1459 FORMAT( ' Metadynamics energy :',E19.10,' (',E15.5,'/ion)') 1653 1460 FORMAT( ' exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1654 1461 FORMAT( ' TAMD energy :',E19.10,' (',E15.5,'/ion)') 1655 1470 FORMAT( ' ion-ion energy :',E19.10,' (',E15.5,'/ion)') 1656 1471 FORMAT(/' Kinetic energy (elc) :',E19.10,' (',E15.5,'/elc)') 1657 1472 FORMAT( ' Kinetic energy (ion) :',E19.10,' (',E15.5,'/ion)') 1658 1473 FORMAT( ' thermostat energy (elc) :',E19.10,' (',E15.5,'/elc)') 1659 1474 FORMAT( ' thermostat energy (ion) :',E19.10,' (',E15.5,'/ion)') 1660 1480 FORMAT(' Temperature : ',F10.1,' K (ion)') 1661 1490 FORMAT(' : ',F10.1,' K (c.o.m.)') 1662 1491 FORMAT(' Temperature : ',F10.1,' K (elc)') 1663 1492 FORMAT(/' Vaverage Eaverage : ',E19.10, E19.10) 1664 1493 FORMAT( ' Vvariance Evariance: ',E19.10, E19.10) 1665 1494 FORMAT( ' Cv - f*kb/(2*nion) : ',E19.10) 1666 1499 FORMAT( ' K.S. SIC-hartree energy :',E19.10, 1667 > ' (',E15.5,'/electron)') 1668 1501 FORMAT( ' K.S. SIC-exc-corr energy :',E19.10, 1669 > ' (',E15.5,'/electron)') 1670 1671 1500 FORMAT(/' orbital energies:') 1672 1505 FORMAT( ' APC energy :',E19.10,' (',E15.5,'/ion)') 1673 1510 FORMAT(2(E18.7,' (',F8.3,'eV)')) 1674 1511 FORMAT(2(E18.7,' (',F8.3,'eV) occ=',F6.3)) 1675 1600 FORMAT(/' Total PSPW energy :',E19.10) 1676 1677 1700 FORMAT(/' QM/MM-pol-vib/CAV Energies') 1678 1701 FORMAT( ' --------------------------') 1679 1702 FORMAT( ' LJ energy :',E19.10) 1680 1703 FORMAT( ' Residual Coulomb energy:',E19.10) 1681 1704 FORMAT( ' MM Vibration energy :',E19.10) 1682 1705 FORMAT( ' QM/MM coupling energy :',E19.10, 1683 > ' (average=',E19.10,' variance=',E19.10,')'/) 1684 1706 FORMAT( ' Average and Variance of QM/MM coupling energy :', 1685 > E19.10,E19.10) 1686 1707 FORMAT( ' QM/MM coupling param. :',E19.10) 1687 1688 1720 FORMAT(/' Dispersion energy :',E19.10) 1689 1690 1691 1740 FORMAT(/' extended Born solvation energies:') 1692 1741 FORMAT(5x,' screen=(epsilon-1)/(epsilon):',F11.6) 1693 1745 FORMAT(5x,' solvation energy (w/o QM polarization) :',E19.10, 1694 > ' (',F8.3,' kcal/mol)') 1695 1746 FORMAT(5x,' solvation energy (w/ QM polarization) :',E19.10, 1696 > ' (',F8.3,' kcal/mol)') 1697 1698 1800 FORMAT(/' Charge Field Energies') 1699 1801 FORMAT( ' ---------------------') 1700 1802 FORMAT( ' - Charge Field/Electron :',E19.10) 1701 1803 FORMAT( ' - Charge Field/Ion :',E19.10) 1702 1804 FORMAT( ' - Charge Field/Charge Field:',E19.10) 1703 1805 FORMAT( ' Charge Field Energy :',E19.10) 1704 1705 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 1706 1707 9000 if (taskid.eq.MASTER) write(luout,9010) ierr 1708 call Parallel2d_Finalize() 1709 call Parallel_Finalize() 1710 1711 cpmdv5 = value 1712 return 1713 END 1714 1715c 1716c Now in nwpw/utilities 1717c 1718c subroutine print_elapsed_time(autime) 1719c implicit none 1720c real*8 autime 1721c 1722c#include "stdio.fh" 1723c 1724c real*8 sectime 1725c 1726c sectime = autime*2.41889d-17 1727c 1728c if (sectime.lt.1.0d-12) then 1729c write(luout,1800) (sectime/1.0d-15)," fs" 1730c else if (sectime.lt.1.0d-9) then 1731c write(luout,1800) (sectime/1.0d-12)," ps" 1732c else 1733c write(luout,1800) (sectime/1.0d-9 )," ns" 1734c end if 1735c 1736c return 1737c 1800 format(//' Elapsed time of simulation was',F8.3,A) 1738c end 1739 1740 1741