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