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