1* 2* $Id$ 3* 4 5*********************************************************************** 6* cgmdv5 * 7 8* This is a developing Born-Oppenheimer MD code for NWCHEM * 9* * 10*********************************************************************** 11 12 logical function cgmdv5(rtdb,flag) 13 implicit none 14 integer rtdb 15 integer flag 16 17#include "global.fh" 18#include "errquit.fh" 19#include "bafdecls.fh" 20#include "btdb.fh" 21#include "inp.fh" 22#include "util.fh" 23#include "stdio.fh" 24 25#include "nwpw_timing.fh" 26 27* *** local variables and parameters **** 28 double precision kb 29 parameter (kb=3.16679d-6) 30 real*8 autoatm 31 parameter (autoatm =290.360032539d6) 32 33 34* **** parallel variables **** 35 integer taskid,np,np_i,np_j 36 integer MASTER 37 parameter(MASTER=0) 38 39* **** timing variables **** 40 real*8 cpu1,cpu2,cpu3,cpu4 41 real*8 t1,t2,t3,t4,av 42 43* **** lattice variables **** 44 integer ngrid(3),nwave,nfft3d,n2ft3d 45 real*8 a,b,c,alpha,beta,gamma 46 47 48* ***** energy variables **** 49 integer ne(2),ms,ispin 50 real*8 E(40),en(2) 51 real*8 dipole(3),cdipole(3) 52 real*8 stress(3,3),lstress(6) 53 54* **** gradient variables **** 55 integer fion(2),fion1(2) 56 57* **** error variables **** 58 logical value,newpsi 59 integer ierr 60 61* **** local variables **** 62 logical verlet,mulliken,SA,calc_pressure,nose,bo_cpmd,vverlet 63 logical oprint,lprint,hprint,qmmm,found,notfirststep,found_bak 64 real*8 gx,gy,gz,cx,cy,cz 65 real*8 vgx,vgy,vgz,vcx,vcy,vcz,ekg,eki0,eki1,dt,fmass,dte 66 real*8 sum,w,emotion_time_shift,wa,wr 67 real*8 EV,EV0,pi,esum1,esum2,eave,evar,dr,E25,E26 68 real*8 ratio,aratio 69 real*8 icharge,cv 70 real*8 sa_alpha(2),sa_decay(2),tt1,tt2,ssr,r,p1,p2,pressure 71 real*8 Te_new,Tr_new,Te_init,Tr_init 72 integer i,k,ia,nion,minimizer,mapping,icount_shift,icount 73 integer nbq,mapping1d,it_in,it_out 74 75 character*50 filename 76 character*255 full_filename,full_bak 77 78 79* **** external functions **** 80 real*8 psp_zv,psp_rc,ewald_rcut 81 real*8 ewald_mandelung 82 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 83 real*8 lattice_unitg,ion_amass 84 integer ewald_ncut,ewald_nshl3d 85 integer psp_lmmax,psp_lmax,psp_locp,ion_nkatm 86 integer psp_nprj,psp_psp_type 87 character*4 ion_atom,ion_aname 88 external psp_zv,psp_rc,ewald_rcut 89 external ewald_mandelung 90 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 91 external lattice_unitg,ion_amass 92 external ewald_ncut,ewald_nshl3d 93 external psp_lmmax,psp_lmax,psp_locp,ion_nkatm 94 external psp_nprj,psp_psp_type 95 external ion_atom,ion_aname 96 97 external psp_comment 98 character*255 psp_comment,comment 99 100 real*8 control_tole,control_tolc,control_tolr,ion_rion,ion_vion 101 external control_tole,control_tolc,control_tolr,ion_rion,ion_vion 102 real*8 control_time_step,control_fake_mass,control_bo_fake_mass 103 external control_time_step,control_fake_mass,control_bo_fake_mass 104 logical control_read,control_move,ion_init,control_out_of_time 105 external control_read,control_move,ion_init,control_out_of_time 106 logical control_translation,control_rotation,control_parallel_io 107 external control_translation,control_rotation,control_parallel_io 108 integer pack_nwave_all 109 integer control_it_in,control_it_out,control_gga,control_version 110 integer control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm 111 integer pspw_charge_nion 112 external pack_nwave_all 113 external control_it_in,control_it_out,control_gga,control_version 114 external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm 115 external pspw_charge_nion 116 117 character*12 control_boundry 118 external control_boundry 119 character*50 control_cell_name 120 external control_cell_name 121 122 logical psp_semicore 123 real*8 psp_rcore,psp_ncore,psp_rlocal 124 external psp_semicore 125 external psp_rcore,psp_ncore,psp_rlocal 126 127 logical psi_initialize, psi_finalize 128 integer psi_ispin,psi_ne,electron_count,linesearch_count 129 external psi_initialize, psi_finalize 130 external psi_ispin,psi_ne,electron_count,linesearch_count 131 real*8 cgsd_energy,cgsd_noit_energy 132 external cgsd_energy,cgsd_noit_energy 133 logical control_Mulliken,control_DOS 134 external control_Mulliken,control_DOS 135 real*8 ion_TotalCharge 136 external ion_TotalCharge 137 logical control_check_charge_multiplicity 138 external control_check_charge_multiplicity 139 logical pspw_charge_found,ion_q_FixIon,ion_q_xyzFixIon 140 external pspw_charge_found,ion_q_FixIon,ion_q_xyzFixIon 141 integer control_minimizer,control_scf_algorithm 142 external control_minimizer,control_scf_algorithm 143 integer control_ks_algorithm 144 external control_ks_algorithm 145 real*8 control_ks_alpha,control_kerker_g0 146 external control_ks_alpha,control_kerker_g0 147 logical control_print,control_balance 148 external control_print,control_balance 149 integer control_mapping,control_mapping1d,control_np_orbital 150 external control_mapping,control_mapping1d,control_np_orbital 151 integer control_ks_maxit_orb,control_ks_maxit_orbs 152 external control_ks_maxit_orb,control_ks_maxit_orbs 153 character*14 ion_q_xyzFixIon_label 154 external ion_q_xyzFixIon_label 155 156 logical pspw_bqext,control_bo_cpmd 157 external pspw_bqext,control_bo_cpmd 158 159 integer control_bo_steps_in,control_bo_steps_out 160 integer control_bo_algorithm 161 integer ion_scaling_atoms_natms,ion_scaling_atoms 162 real*8 control_bo_time_step,control_rti,ion_ke,ion_Temperature 163 external control_bo_steps_in,control_bo_steps_out 164 external control_bo_algorithm 165 external ion_scaling_atoms_natms,ion_scaling_atoms 166 external control_bo_time_step,control_rti,ion_ke,ion_Temperature 167 real*8 ion_com_Temperature 168 external ion_com_Temperature 169 logical control_Nose,control_SA,control_pressure,Nose_restart 170 real*8 control_Nose_Te,control_Nose_Tr,control_SA_decay,Nose_ssr 171 real*8 Nose_Pr,Nose_Pe,Nose_Ee0,Nose_Er0,Nose_dXr 172 real*8 Nose_Tr,Nose_Te,Nose_Qr,Nose_Qe,Nose_r_energy 173 integer Nose_Mchain,Nose_Nchain 174 external control_Nose,control_SA,control_pressure,Nose_restart 175 external control_Nose_Te,control_Nose_Tr,control_SA_decay,Nose_ssr 176 external Nose_Pr,Nose_Pe,Nose_Ee0,Nose_Er0,Nose_dXr 177 external Nose_Tr,Nose_Te,Nose_Qr,Nose_Qe,Nose_r_energy 178 external Nose_Mchain,Nose_Nchain 179 integer ion_nconstraints,ion_ndof 180 external ion_nconstraints,ion_ndof 181 logical ion_makehmass2,control_periodic_dipole 182 external ion_makehmass2,control_periodic_dipole 183 logical control_precondition,control_fractional 184 external control_precondition,control_fractional 185 186 187* |************| 188*****************************| PROLOGUE |**************************** 189* |************| 190 191 value = .true. 192 pi = 4.0d0*datan(1.0d0) 193 194 call nwpw_timing_init() 195 call dcopy(30,0.0d0,0,E,1) 196 197 198* **** get parallel variables **** 199 call Parallel_Init() 200 call Parallel_np(np) 201 call Parallel_taskid(taskid) 202 203 value = control_read(11,rtdb) 204 if (.not. value) 205 > call errquit('error reading control',0, INPUT_ERR) 206 207 call Parallel2d_Init(control_np_orbital()) 208 call Parallel2d_np_i(np_i) 209 call Parallel2d_np_j(np_j) 210 211 oprint= ((taskid.eq.MASTER).and.control_print(print_medium)) 212 lprint= ((taskid.eq.MASTER).and.control_print(print_low)) 213 hprint= ((taskid.eq.MASTER).and.control_print(print_high)) 214 215 216 if (oprint) call current_second(cpu1) 217 218* ***** print out header **** 219 if (oprint) then 220 write(luout,1000) 221 write(luout,1010) 222 write(luout,1020) 223 write(luout,1010) 224 write(luout,1030) 225 write(luout,1031) 226 write(luout,1010) 227 write(luout,1035) 228 write(luout,1010) 229 write(luout,1040) 230 write(luout,1010) 231 write(luout,1041) 232 write(luout,1010) 233 write(luout,1000) 234 call nwpw_message(1) 235 write(luout,1110) 236 end if 237 238 ngrid(1) = control_ngrid(1) 239 ngrid(2) = control_ngrid(2) 240 ngrid(3) = control_ngrid(3) 241 nwave = 0 242 minimizer = control_minimizer() 243 mapping = control_mapping() 244 245* **** initialize pspw_director **** 246 call pspw_director_init(rtdb) 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 257* ***** Initialize double D3dB data structure **** 258 if (control_version().eq.4) 259 > call D3dB_Init(2,2*ngrid(1),2*ngrid(2),2*ngrid(3),mapping) 260 261 262* **** initialize lattice data structure **** 263 call lattice_init() 264 call G_init() 265 call mask_init() 266 call Pack_init() 267 call D3dB_pfft_init() 268 call ga_sync() 269 270* **** read ions **** 271 value = ion_init(rtdb) 272 273* **** initialize FixIon constraint **** 274 call ion_init_FixIon(rtdb) 275 276* **** allocate psp data structure and read in psedupotentials into it **** 277 call psp_init() 278 call psp_readall() 279 if (psp_semicore(0)) call semicore_check() 280 281* **** initialize G,mask,ke,coulomb,and ewald data structures **** 282 call ke_init() 283 if (control_version().eq.3) call coulomb_init() 284 if (control_version().eq.4) call coulomb2_init() 285 call strfac_init() 286 call phafac() 287 if (control_version().eq.3) then 288 call ewald_init() 289 call ewald_phafac() 290 end if 291 292* **** read in wavefunctions and initialize psi **** 293 if (.not.control_check_charge_multiplicity()) then 294 call psi_new() 295 newpsi = .true. 296 else 297 newpsi = .false. 298 end if 299 300* **** Initialize 1d mappings for ne spaces **** 301 call psi_get_ne(ispin,ne) 302 mapping1d = control_mapping1d() 303 call Dne_init(ispin,ne,mapping1d) 304 305 306* **** read in wavefunctions and initialize psi **** 307 value = psi_initialize() 308 309 310 !call D3dB_n_fft_init(1,psi_ne(1)+psi_ne(2)) 311 312* **** electron and geodesic data structures **** 313 call electron_init() 314 if ((minimizer.eq.4).or.(minimizer.eq.7)) then !debug 315 call geodesic2_init() 316 else 317 call geodesic_init() 318 end if 319 320* **** initialize QM/MM **** 321 call pspw_init_APC(rtdb) 322 call pspw_qmmm_init(rtdb) 323 call pspw_charge_init(rtdb) 324 325* **** initialize SIC and HFX **** 326 ne(1) = psi_ne(1) 327 ne(2) = psi_ne(2) 328 call pspw_init_SIC(rtdb,ne) 329 call pspw_init_HFX(rtdb,psi_ispin(),ne) 330 331* **** initialize DFT+U **** 332 call psp_U_init() 333 334* **** initialize Meta GGA **** 335 call nwpw_meta_gga_init(control_gga()) 336 337* **** initialize metadynamics **** 338 call meta_initialize(rtdb) 339 340 341 342 343* **** initialize pressure **** 344 calc_pressure = control_pressure().and.(control_version().eq.3) 345 pressure = 0.0d0 346 p1 = 0.0d0 347 p2 = 0.0d0 348 if (calc_pressure) then 349 call psp_stress_init() 350 call psp_stress_readall() 351 end if 352 353 354* **** initialize frac_occ data structure **** 355c call frac_occ_init(rtdb,psi_ispin(),ne) 356 357* **** initialize linesearching **** 358 call linesearch_init() 359 360 361* ****************************** 362* **** scaling ion velocity **** 363* ****************************** 364 call ion_init_ke(ekg,eki0,eki1) 365 366 367* **** Initialize thermostats **** 368 nose = control_Nose() 369 if (nose) call Nose_Init((ne(1)+ne(2)),0.01d0) 370 371* **** Initialize simulated annealing **** 372 SA = .false. 373 Tr_init = 0.0d0 374 sa_alpha(2) = 1.0d0 375 if (control_SA()) then 376 if (nose) then 377 SA = .true. 378 sa_decay(2) = control_SA_decay(2) 379 Tr_init = control_Nose_Tr() 380 else 381 dt = control_bo_time_step() 382 SA = .false. 383 sa_decay(2) = control_SA_decay(2) 384 sa_alpha(2) = dexp( -(dt/control_SA_decay(2)) ) 385 end if 386 end if 387 388 bo_cpmd = control_bo_cpmd() 389 vverlet = (control_bo_algorithm().eq.1) 390 391 392* |**************************| 393****************** summary of input data ********************** 394* |**************************| 395 call center_geom(cx,cy,cz) 396 call center_mass(gx,gy,gz) 397 call center_v_geom(vcx,vcy,vcz) 398 call center_v_mass(vgx,vgy,vgz) 399 mulliken = control_Mulliken() 400 401 402 if (oprint) then 403 write(luout,1111) np 404 write(luout,1117) np_i,np_j 405 if (mapping.eq.1) write(luout,1112) 406 if (mapping.eq.2) write(luout,1113) 407 if (mapping.eq.3) write(luout,1118) 408 if (control_balance()) then 409 write(luout,1114) 410 else 411 write(luout,1116) 412 end if 413 write(luout,1115) 414 if (control_parallel_io()) then 415 write(luout,1123) 416 else 417 write(luout,1124) 418 end if 419 420 write(luout,1121) control_boundry(),control_version() 421 if (psi_ispin().eq.1) write(luout,1130) "restricted" 422 if (psi_ispin().eq.2) write(luout,1130) "unrestricted" 423 !if (qmmm) write(luout,1122) 424 if (control_fractional()) write(luout,1132) 425 426 427 call v_bwexc_print(luout,control_gga()) 428 429 call pspw_print_SIC(luout) 430 call pspw_print_HFX(luout) 431 call nwpw_meta_gga_print(luout) 432 if (ion_makehmass2()) write(luout,1135) 433 write(luout,1140) 434 do ia = 1,ion_nkatm() 435 call psp_print(ia) 436c call psp_check_print(ia) 437 end do 438 439 440 441 icharge = -(psi_ne(1)+psi_ne(psi_ispin())) 442 en(1) = psi_ne(1) 443 en(psi_ispin()) = psi_ne(psi_ispin()) 444 445c if (fractional) then 446c icharge = 0 447c do ms=1,psi_ispin() 448c en(ms) =0.0 449c do i=1,ne(ms) 450c k = fweight(1)+(i-1)+(ms-1)*ne(1) 451c icharge = icharge - (3-psi_ispin())*dbl_mb(k) 452c en(ms) = en(ms) + dbl_mb(k) 453c end do 454c end do 455c end if 456 icharge = icharge + ion_TotalCharge() 457 write(luout,1159) icharge 458 459 write(luout,1160) 460 write(luout,1170) (ion_atom(K),ion_natm(K),K=1,ion_nkatm()) 461 !if (hprint) then 462 write(luout,1180) 463 do I=1,ion_nion() 464 if (ion_q_FixIon(I)) then 465 write(luout,1191) I,ion_aname(I),(ion_rion(K,I),K=1,3), 466 > ion_amass(I)/1822.89d0 467 else if (ion_q_xyzFixIon(I)) then 468 write(luout,1194) I,ion_aname(I),(ion_rion(K,I),K=1,3), 469 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 470 else 471 write(luout,1190) I,ion_aname(I),(ion_rion(K,I),K=1,3), 472 > ion_amass(I)/1822.89d0 473 end if 474 end do 475 write(luout,1200) cx,cy,cz 476 write(luout,1210) gx,gy,gz 477 478 call pspw_charge_Print(luout) 479 480 write(luout,1181) 481 write(luout,1192) (I,ion_aname(I), 482 > (ion_vion(K,I),K=1,3),I=1,ion_nion()) 483 write(luout,1200) vcx,vcy,vcz 484 write(luout,1210) vgx,vgy,vgz 485 write(luout,1211) ion_nconstraints(),ion_ndof() 486 487 !end if 488 489 !write(6,1220) psi_ne(1),psi_ne(psi_ispin()),' (Fourier space)' 490c if (fractional) then 491c write(6,1219) en(1),en(psi_ispin()),' ( fractional)' 492c write(6,1221) psi_ne(1),psi_ne(psi_ispin()),' (Fourier space)' 493c else 494 write(luout,1220) psi_ne(1),psi_ne(psi_ispin()), 495 > ' (Fourier space)' 496 write(luout,1221) psi_ne(1),psi_ne(psi_ispin()), 497 > ' (Fourier space)' 498c end if 499 write(luout,1230) 500 write(luout,1233) control_cell_name() 501 write(luout,1241) lattice_unita(1,1), 502 > lattice_unita(2,1), 503 > lattice_unita(3,1) 504 write(luout,1242) lattice_unita(1,2), 505 > lattice_unita(2,2), 506 > lattice_unita(3,2) 507 write(luout,1243) lattice_unita(1,3), 508 > lattice_unita(2,3), 509 > lattice_unita(3,3) 510 write(luout,1244) lattice_unitg(1,1), 511 > lattice_unitg(2,1), 512 > lattice_unitg(3,1) 513 write(luout,1245) lattice_unitg(1,2), 514 > lattice_unitg(2,2), 515 > lattice_unitg(3,2) 516 write(luout,1246) lattice_unitg(1,3), 517 > lattice_unitg(2,3), 518 > lattice_unitg(3,3) 519 call lattice_abc_abg(a,b,c,alpha,beta,gamma) 520 write(luout,1232) a,b,c,alpha,beta,gamma 521 write(luout,1231) lattice_omega() 522 write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 523 > pack_nwave_all(0),pack_nwave(0) 524 write(luout,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 525 > pack_nwave_all(1),pack_nwave(1) 526 if (control_version().eq.3) then 527 write(luout,1260) ewald_rcut(),ewald_ncut() 528 write(luout,1261) ewald_mandelung() 529 end if 530 531 write(luout,1270) 532 write(luout,1280) control_time_step(),control_fake_mass() 533 write(luout,1290) control_tole(),control_tolc() 534 write(luout,1281) control_it_in()*control_it_out(), 535 > control_it_in(),control_it_out() 536 if ((minimizer.eq.5).or.(minimizer.eq.8)) then 537 write(luout,1291) 538 539 if (control_ks_algorithm().eq.-1) then 540 if (control_precondition()) then 541 write(luout,1292) 542 > "block preconditioned conjugate gradient" 543 else 544 write(luout,1292) "block conjugate gradient" 545 end if 546 end if 547 if (control_ks_algorithm().eq.0) then 548 if (control_precondition()) then 549 write(luout,1292) "preconditioned conjugate gradient" 550 else 551 write(luout,1292) "conjugate gradient" 552 end if 553 end if 554 if (control_ks_algorithm().eq.1) 555 > write(luout,1292) "rmm-diis" 556 write(luout,1295) control_ks_maxit_orb(), 557 > control_ks_maxit_orbs() 558 if (control_scf_algorithm().eq.0) 559 > write(luout,1293) "simple mixing" 560 if (control_scf_algorithm().eq.1) 561 > write(luout,1293) "Anderson potential mixing" 562 if (control_scf_algorithm().eq.2) 563 > write(luout,1293) "Johnson-Pulay mixing" 564 if (control_scf_algorithm().eq.3) 565 > write(luout,1293) "Anderson density mixing" 566 write(luout,1294) control_ks_alpha() 567 write(luout,1301) control_kerker_g0() 568 end if 569 write(luout,1310) 570 if (.not.control_translation()) write(luout,1271) 571 if (.not.control_rotation()) write(luout,1272) 572 if (.not.bo_cpmd) then 573 write(luout,1320) control_bo_time_step(), 574 > control_bo_steps_in()*control_bo_steps_out(), 575 > control_bo_steps_in(),control_bo_steps_out() 576 else 577 write(luout,1321) control_bo_time_step(), 578 > control_bo_fake_mass(), 579 > control_bo_steps_in()*control_bo_steps_out(), 580 > control_bo_steps_in(),control_bo_steps_out() 581 end if 582 583 if(control_bo_algorithm().eq.0) write(6,1330) "Position Verlet" 584 if(control_bo_algorithm().eq.1) write(6,1330) "Velocity Verlet" 585 if(control_bo_algorithm().eq.2) write(6,1330) "Leap Frog" 586 587 write(luout,1340) control_rti() 588 call ion_scaling_atoms_print(luout) 589 write(luout,1222) eki0,ekg 590 write(luout,1223) eki1 591 write(luout,1224) (eki1-eki0) 592 if (nose) then 593 write(luout,1395) 594 if (Nose_restart()) then 595 write(luout,*) " thermostats resused" 596 else 597 write(luout,*) " thermostats initialized" 598 end if 599 do i=1,Nose_Nchain() 600 write(luout,1398) i,control_Nose_Tr(),Nose_Qr(i), 601 > Nose_Pr(i),Nose_Er0(i) 602 end do 603 else 604 write(luout,1394) 605 end if 606 if (calc_pressure) write(luout,1393) 607 if (control_SA()) then 608 write(luout,1396) sa_decay(2) 609 end if 610 if (mulliken) write(luout,1399) 611 write(luout,1300) 612 call util_flush(luout) 613 end if 614 615* |***************************| 616****************** simple MD loop ********************** 617* |***************************| 618 if (taskid.eq.MASTER) call current_second(cpu2) 619 if ((taskid.eq.MASTER).and.(.not.calc_pressure)) 620 > call nwpw_message(10) 621 if ((taskid.eq.MASTER).and.(calc_pressure)) 622 > call nwpw_message(11) 623 624 625* **** write initial position to xyz data **** 626 call xyz_init() ! unit=18 627 call MOTION_init(rtdb) ! unit=19 628 629 630* *** fei io **** 631 call fei_init(rtdb) 632 633* **** ecce print **** 634 call ecce_print_module_entry('task Born-Oppenheimer') 635 !call ecce_print_module_entry('driver') 636 call movecs_ecce_print_off() 637 638 639* **** finalize pressure **** 640 if (calc_pressure) then 641 call psp_stress_end() 642 end if 643 644 645 646 647* ************************************ 648* **** open up other MOTION files **** 649* ************************************ 650 651* **** open EMOTION file **** 652 E25 = 0.0d0 653 E26 = 0.0d0 654 icount_shift = 0 655 if (.not.btdb_cget(rtdb,'nwpw:emotion_filename',1,filename)) 656 > call util_file_prefix('emotion',filename) 657 call util_file_name_noprefix(filename,.false., 658 > .false., 659 > full_filename) 660 if (taskid.eq.MASTER) then 661 662* **** check for backup file **** 663 call util_file_name_noprefix('EMOTION99-bak',.false., 664 > .false., 665 > full_bak) 666 inquire(file=full_bak,exist=found_bak) 667 if (found_bak) then 668 write(*,*) 669 write(*,*) "EMOTION99-bak exists:" 670 i=index(full_bak,' ') 671 k=index(full_filename,' ') 672 write(*,*) " Copying ",full_bak(1:i), 673 > " to ",full_filename(1:k) 674 write(*,*) 675 call util_file_copy(full_bak,full_filename) 676 end if 677 678 emotion_time_shift = 0.0d0 679 icount_shift = 0 680 inquire(file=full_filename,exist=found) 681 if (found) then 682 open(unit=31,file=full_filename,form='formatted', 683 > status='old') 684 do while (found) 685 read(31,*,end=100) emotion_time_shift,w,sum 686 E25 = E25 + sum !*** take care of running sums *** 687 E26 = E26 + sum*sum 688 icount_shift = icount_shift + 1 689 end do 690 100 continue 691#if defined(FUJITSU_SOLARIS) || defined(PSCALE) || defined(__crayx1) || defined(GCC46) 692 backspace 31 693#endif 694 else 695 open(unit=31,file=full_filename,form='formatted', 696 > status='new') 697 end if 698 end if 699 700* **** open EIGMOTION file **** 701 if (mulliken) then 702 if (.not.btdb_cget(rtdb,'nwpw:eigmotion_filename',1,filename)) 703 > call util_file_prefix('eigmotion',filename) 704 call util_file_name_noprefix(filename,.false., 705 > .false., 706 > full_filename) 707 if (taskid.eq.MASTER) 708 > open(unit=32,file=full_filename,form='formatted') 709 end if 710 711 call xyz_write() 712 713* **** allocate fion *** 714 nion = ion_nion() 715 if (pspw_charge_found().and. 716 > (.not.pspw_bqext())) nion = nion + pspw_charge_nion() 717 718 value = BA_push_get(mt_dbl,(3*nion), 719 > 'fion',fion(2),fion(1)) 720 if (vverlet) then 721 value = value.and.BA_push_get(mt_dbl,(3*nion), 722 > 'fion1',fion1(2),fion1(1)) 723 end if 724 if (.not. value) call errquit('out of stack memory',0, MA_ERR) 725 726 727 if (bo_cpmd) call psi_cpmd_start() 728 729* ***** do a simple md loop *** 730 notfirststep = .false. 731 dt = control_bo_time_step() 732 fmass = control_bo_fake_mass() 733 dte = dt*dt/fmass 734 call control_reduce_print() 735 EV = cgsd_energy(newpsi) 736 call cgsd_energy_gradient_md(dbl_mb(fion(1))) 737 738 r = 1.0d0 739 if (nose) r = (1.0d0-0.5d0*dt*Nose_dXr()) 740 call ion_newton_step(dbl_mb(fion(1)),sa_alpha(2)*r) 741 eki1 = ion_ke() 742 if (nose) call Nose_Newton_Step(0.01d0,eki1) 743 744 call xyz_write() 745 746 it_out = control_bo_steps_out() 747 it_in = control_bo_steps_in() 748 icount = 0 749 if (it_out.lt.1) goto 102 750 751 752 101 continue 753 icount = icount + 1 754 755c **** inner loop **** 756 do i=1,it_in 757 if (vverlet) then 758 call Parallel_shared_vector_copy(.true.,3*nion, 759 > dbl_mb(fion(1)), 760 > dbl_mb(fion1(1))) 761 call ion_shift21() 762 else 763 call ion_shift() 764 if (nose) call Nose_shift() 765 end if 766 767 if (bo_cpmd) then 768 if (notfirststep) then 769 call psi_cpmd_step(dte) 770 else 771 notfirststep = .true. 772 call psi_cpmd_step(0.5d0*dte) 773 end if 774 end if 775 776 EV = cgsd_energy(.false.) 777 call cgsd_energy_gradient_md(dbl_mb(fion(1))) 778 779 780 if (nose) then 781 ssr = Nose_ssr() 782 call ion_nose_step(ssr,dbl_mb(fion(1))) 783 eki1 = ion_ke() 784 call Nose_Verlet_Step(0.01d0,eki1) 785 else if (vverlet) then 786 call ion_vverlet_step(dbl_mb(fion(1)), 787 > dbl_mb(fion1(1))) 788 call ion_vshift() 789 call ion_newton_step(dbl_mb(fion(1)),sa_alpha(2)) 790 eki1 = ion_ke() 791 else 792 call ion_verlet_step(dbl_mb(fion(1)),sa_alpha(2)) 793 eki1 = ion_ke() 794 end if 795 end do 796 E25 = E25 + EV 797 E26 = E26 + EV*EV 798 799 800 if (calc_pressure) then 801 call psi_1pressure_stress(pressure,p1,p2,stress) 802 end if 803 804 if (oprint) then 805 if (calc_pressure) then 806 if (nose) then 807 write(luout,1350) icount*it_in,EV+eki1+Nose_r_energy(), 808 > EV,eki1,ion_Temperature(), 809 > pressure*autoatm 810 write(31,1360) icount*it_in*dt+emotion_time_shift, 811 > EV+eki1+Nose_r_energy(), 812 > EV,eki1,ion_Temperature(), 813 > pressure*autoatm 814 else 815 write(luout,1350) icount*it_in,EV+eki1,EV,eki1, 816 > ion_Temperature(), 817 > pressure*autoatm 818 write(31,1360) icount*it_in*dt+emotion_time_shift, 819 > EV+eki1,EV,eki1,ion_Temperature(), 820 > pressure*autoatm 821 end if 822 else 823 if (nose) then 824 write(luout,1350) icount*it_in,EV+eki1+Nose_r_energy(), 825 > EV,eki1,ion_Temperature() 826 write(31,1360) icount*it_in*dt+emotion_time_shift, 827 > EV+eki1+Nose_r_energy(), 828 > EV,eki1,ion_Temperature() 829 else 830 write(luout,1350) icount*it_in,EV+eki1,EV,eki1, 831 > ion_Temperature() 832 write(31,1360) icount*it_in*dt+emotion_time_shift, 833 > EV+eki1,EV,eki1,ion_Temperature() 834 end if 835 end if 836 call util_flush(6) 837 call util_flush(31) 838 end if 839 840 call fei_output(EV,dbl_mb(fion(1))) 841 call xyz_write() 842 call MOTION_write(icount*it_in*dt) 843 844* **** call pspw_director **** 845 call pspw_director(rtdb) 846 847* **** update thermostats using SA decay **** 848 if (SA) then 849 tt2 = icount*it_in*dt/sa_decay(2) 850 Tr_new = Tr_init*dexp(-tt2) 851 call Nose_reset_T(Tr_new,Tr_new) 852 end if 853 854* **** exit early **** 855 if (control_out_of_time()) then 856 if (oprint) 857 > write(luout,*) ' *** out of time. iteration terminated' 858 go to 102 859 end if 860 861 862 if (icount.lt.it_out) go to 101 863 if (oprint) write(luout,*) 864 > '*** arrived at the Maximum iteration. terminated.' 865 866*:::::::::::::::::::: end of iteration loop ::::::::::::::::::::::::: 867 call control_up_print() 868 869 870 102 continue 871 872 if (bo_cpmd) call psi_cpmd_end() 873 874* **** close xyz and MOTION files **** 875 call xyz_end() 876 call MOTION_end() 877 if (taskid.eq.MASTER) then 878 close(unit=31) 879 if (mulliken) close(unit=32) 880 881* **** remove EMOTION backup file *** 882 call util_file_name_noprefix('EMOTION99-bak',.false., 883 > .false., 884 > full_bak) 885 call util_file_unlink(full_bak) 886 887 end if 888 889* *** close fei io **** 890 call fei_end() 891 892 893 894c* ***** do a simple mc loop *** 895c call control_reduce_print() 896c dr = 0.08d0 897c aratio = 0.5d0 898c wa = 0.0d0 899c wr = 0.0d0 900c esum1 = 0.0d0 901c esum2 = 0.0d0 902c eave = 0.0d0 903c evar = 0.0d0 904c beta = util_random(39) !**seed the random number generator 905c beta = 1.0d0/(kb*298.15) 906c call xyz_write() 907c EV0 = cgsd_energy(newpsi) 908c do it=1,control_bo_steps() 909c 910c call ion_MC_step(dr) 911c EV = cgsd_energy(.false.) 912c w = dexp(-beta*(EV-EV0)) 913c alpha = util_random(0) 914c 915c !*** accept the sequence *** 916c if (alpha.lt.w) then 917c EV0 = EV 918c !write(6,*) " ----ACCEPTED----" 919c if (taskid.eq.MASTER) then 920c write(6,1350) it,EV,EV0,w,alpha,wa,wr,eave,evar,dr 921c call util_flush(6) 922c end if 923c call xyz_write() 924c call MOTION_write(it*dt) 925c wa = wa + 1.0d0 926c esum1 = esum1 + EV 927c esum2 = esum2 + EV*EV 928c eave = esum1/wa 929c evar = esum2/wa 930c evar = evar - eave*eave 931c 932c !*** reject the sequence *** 933c else 934c call ion_MC_reject_step() 935c !write(6,*) " ----REJECTED----" 936c wr = wr + 1.0d0 937c end if 938c 939c !*** adjust *** 940c if (mod(it,10).eq.0) then 941c ratio = wa/(wa+wr) 942c dr = dr*(ratio/aratio) 943c !if (ratio.lt.aratio) dr = dr*0.95d0 944c !if (ratio.gt.aratio) dr = dr*1.05d0 945c end if 946c 947c end do 948c 949c call control_up_print() 950 951* |***************************| 952****************** simple MD loop ********************** 953* |***************************| 954 955 956* |***************************| 957****************** report summary of results ********************** 958* |***************************| 959 960 call center_geom(cx,cy,cz) 961 call center_mass(gx,gy,gz) 962 call center_v_geom(vcx,vcy,vcz) 963 call center_v_mass(vgx,vgy,vgz) 964 965 if (oprint) then 966 call print_elapsed_time(icount*it_in*dt) 967 write(luout,1300) 968 write(luout,1410) 969 write(luout,1420) 970 do I=1,ion_nion() 971 if (ion_q_FixIon(I)) then 972 write(luout,1191) I,ion_aname(I),(ion_rion(k,i),K=1,3), 973 > ion_amass(I)/1822.89d0 974 else if (ion_q_xyzFixIon(I)) then 975 write(luout,1194) I,ion_aname(I),(ion_rion(k,i),K=1,3), 976 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 977 else 978 write(luout,1190) I,ion_aname(I),(ion_rion(k,i),K=1,3), 979 > ion_amass(I)/1822.89d0 980 end if 981 end do 982 write(luout,1200) cx,cy,cz 983 write(luout,1210) gx,gy,gz 984 985 write(luout,1421) 986 write(luout,1192) (I,ion_aname(I), 987 > (ion_vion(K,I),K=1,3),I=1,ion_nion()) 988 write(luout,1200) vcx,vcy,vcz 989 write(luout,1210) vgx,vgy,vgz 990 write(luout,1211) ion_nconstraints(),ion_ndof() 991 end if 992 993 EV = cgsd_noit_energy() 994 995 if (oprint) then 996 997 write(luout,1300) 998 write(luout,1472) ion_ke(),ion_ke()/ion_nion() 999 1000* **** write out Temperatures **** 1001 eki0 = ion_Temperature() 1002 write(luout,1480) eki0 1003 write(luout,1490) ion_com_Temperature() 1004 1005 eave = E25/dble(icount+icount_shift) 1006 evar = E26/dble(icount+icount_shift) 1007 evar = evar - eave*eave 1008 cv = (evar)/(kb*ion_Temperature()**2) 1009 cv = cv/dble(ion_nion()) 1010 write(luout,1492) eave 1011 write(luout,1493) evar 1012 write(luout,1494) cv 1013 1014 end if 1015 call cgsd_energy_gradient(dbl_mb(fion(1))) 1016 1017* **** calculate excited state orbitals **** 1018 call cgsd_excited() 1019 1020* **** calculate oep eigenvalues **** 1021 call cgsd_oep_eigenvalues() 1022 1023* **** extra energy output for QA test **** 1024 if (lprint) write(6,1600) EV 1025 1026* **** calculate the spin contamination **** 1027 if (flag.gt.-1) call psi_spin2(dipole(1)) 1028 1029* **** calculate the dipole *** 1030 dipole(1) = 0.0d0 1031 dipole(2) = 0.0d0 1032 dipole(3) = 0.0d0 1033 cdipole(1) = 0.0d0 1034 cdipole(2) = 0.0d0 1035 cdipole(3) = 0.0d0 1036 if (flag.gt.-1) then 1037 call rho_dipole(dipole) 1038 if (control_periodic_dipole()) call psi1_crystal_dipole(cdipole) 1039 end if 1040 1041 1042* **** calculate the stress tensor **** 1043 call dcopy(9,0.0d0,0,stress,1) 1044 call dcopy(6,0.0d0,0,lstress,1) 1045 if (flag.eq.3) then 1046 call psp_stress_init() 1047 call psp_stress_readall() 1048 call cgsd_energy_stress(stress,lstress) 1049 call psp_stress_end() 1050 end if 1051 1052 1053* ************************************************************* 1054* **** output energy, dipole, and gradient to rtdb for use **** 1055* **** by task_energy and task_gradient **** 1056* ************************************************************* 1057 if (flag.gt.-1) then 1058 value = btdb_put(rtdb,'pspw:energy',mt_dbl,1,EV) 1059 value = value.and. 1060 > btdb_put(rtdb,'pspw:dipole',mt_dbl, 1061 > 3,dipole) 1062 value = value.and. 1063 > btdb_put(rtdb,'pspw:crystal_dipole',mt_dbl, 1064 > 3,cdipole) 1065 end if 1066 if (flag.gt.0) then 1067 value = value.and. 1068 > btdb_put(rtdb,'pspw:gradient',mt_dbl, 1069 > 3*nion,dbl_mb(fion(1))) 1070 end if 1071 if (flag.eq.3) then 1072 value = value.and. 1073 > btdb_put(rtdb,'pspw:stress',mt_dbl, 1074 > 9,stress) 1075 value = value.and. 1076 > btdb_put(rtdb,'pspw:lstress',mt_dbl, 1077 > 6,lstress) 1078 end if 1079 if (vverlet) value = value.and.BA_pop_stack(fion1(2)) 1080 value = value.and.BA_pop_stack(fion(2)) 1081 if (.not. value) call errquit('cgmdv5: error writing rtdb',0, 1082 & RTDB_ERR) 1083 1084 if (taskid.eq.MASTER) call current_second(cpu3) 1085 1086* |***************************| 1087****************** Epilogue ********************** 1088* |***************************| 1089 1090* **** calculate Mulliken Populations *** 1091 if (control_Mulliken()) call psi_Mulliken(rtdb) 1092 1093 1094 1095 1096* **** calculate Mulliken Populations *** 1097 1098* **** write wavefunctions to file and finalize psi **** 1099 if (flag.eq.-1) then 1100 value = psi_finalize(.false.) 1101 else 1102 value = psi_finalize(.true.) 1103 end if 1104 1105* **** write geometry to rtdb **** 1106 call pspw_charge_write(rtdb) 1107 call ion_write(rtdb) 1108 1109 1110* **** deallocate heap memory **** 1111 call electron_finalize() 1112 if ((minimizer.eq.4).or.(minimizer.eq.7)) then 1113 call geodesic2_finalize() 1114 else 1115 call geodesic_finalize() 1116 end if 1117 if (control_version().eq.3) call ewald_end() 1118 call strfac_end() 1119 if (control_version().eq.3) call coulomb_end() 1120 if (control_version().eq.4) call coulomb2_end() 1121 call ke_end() 1122 call mask_end() 1123 call Pack_end() 1124 call G_end() 1125 call psp_U_end() 1126 call nwpw_meta_gga_end() 1127 call pspw_end_SIC() 1128 call pspw_end_HFX() 1129 call pspw_end_APC() 1130 call pspw_charge_end() 1131 call pspw_qmmm_end() 1132c call frac_occ_end() 1133 call meta_finalize(rtdb) 1134 1135 if (nose) call Nose_end() 1136 call ion_end() 1137 call psp_end() 1138 call ion_end_FixIon() 1139 !call D3dB_n_fft_end(1) 1140 call D3dB_pfft_end() 1141 call D3dB_end(1) 1142 if (control_version().eq.4) call D3dB_end(2) 1143 call Dne_end() 1144 call psi_data_end() 1145 1146 1147* **** do anaylysis on MOTION files **** 1148 call cpmd_properties(rtdb) 1149 1150 1151 1152* |***************************| 1153****************** report consumed cputime ********************** 1154* |***************************| 1155 if (lprint) then 1156 CALL current_second(cpu4) 1157 1158 T1=CPU2-CPU1 1159 T2=CPU3-CPU2 1160 T3=CPU4-CPU3 1161 T4=CPU4-CPU1 1162 AV=T2/dble(electron_count()) 1163 write(6,1801) 1164 write(6,1802) 1165 write(6,1803) T1 1166 write(6,1804) T2 1167 write(6,1805) T3 1168 write(6,1806) T4 1169 write(6,1807) AV,electron_count(),linesearch_count() 1170 1171 call nwpw_timing_print_final(oprint,electron_count()) 1172 1173 write(6,*) 1174 CALL nwpw_MESSAGE(4) 1175 end if 1176 1177 1178 1179 call Parallel_Finalize() 1180 cgmdv5 = value 1181 return 1182 1183 1184*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 1185 1000 FORMAT(10X,'****************************************************') 1186 1010 FORMAT(10X,'* *') 1187 1020 FORMAT(10X,'* NWPW PSPW Calculation *') 1188 1030 FORMAT(10X,'* [ Born-Oppenheimer molecular ] *') 1189 1031 FORMAT(10X,'* [ dynamics simulation ] *') 1190 1035 FORMAT(10x,'* [ NorthWest Chemistry implementation ] *') 1191 1040 FORMAT(10X,'* version #1.00 07/18/06 *') 1192 1041 FORMAT(10X,'* This code was developed by Eric J. Bylaska. *') 1193 1100 FORMAT(//) 1194 1110 FORMAT(10X,'================ input data ========================') 1195 1111 FORMAT(/' number of processors used:',I10) 1196 1112 FORMAT( ' parallel mapping : 1d-slab') 1197 1113 FORMAT( ' parallel mapping : 2d-hilbert') 1198 1114 FORMAT( ' parallel mapping : balanced') 1199 1115 FORMAT(/' options:') 1200 1116 FORMAT( ' parallel mapping : not balanced') 1201 1117 FORMAT( ' processor grid :',I4,' x',I4) 1202 1118 FORMAT( ' parallel mapping : 2d-hcurve') 1203 1120 FORMAT(5X,' ionic motion = ',A) 1204 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 1205 1122 FORMAT(5X,' qmmm simulation') 1206 1123 FORMAT( ' parallel io : on') 1207 1124 FORMAT( ' parallel io : off') 1208 1130 FORMAT(5X,' electron spin = ',A) 1209 1131 FORMAT(5X,' exchange-correlation = ',A) 1210 1132 FORMAT(5X,' using fractional occupation') 1211 1135 FORMAT(/' The masses of QM H atoms converted to 2.0 amu. ', 1212 > /' To turn off this default', 1213 > ' set nwpw:makehmass2 .false.') 1214 1140 FORMAT(/' elements involved in the cluster:') 1215 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I3) 1216 1151 FORMAT(5X,' cutoff =',4F8.3) 1217 1152 FORMAT(12X,' highest angular component : ',i3) 1218 1153 FORMAT(12X,' local potential used : ',i3) 1219 1154 FORMAT(12X,' number of non-local projections: ',i3) 1220 1155 FORMAT(12X,' semicore corrections included : ', 1221 > F6.3,' (radius) ',F6.3,' (charge)') 1222 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 1223 1157 FORMAT(12X,' comment : ',A) 1224 1158 FORMAT(12X,' pseudpotential type : ',i3) 1225 1226 1159 FORMAT(/' total charge:',F8.3) 1227 1160 FORMAT(/' atomic composition:') 1228 1170 FORMAT(7(5X,A4,':',I5)) 1229 1180 FORMAT(/' initial position of ions (au):') 1230 1181 FORMAT(/' initial velocity of ions after scaling (au):') 1231 1232 1190 FORMAT(5X, I4, A5 ,' (',3F11.5,' ) - atomic mass= ',F7.3,' ') 1233 1191 FORMAT(5X, I4, A5, ' (',3F11.5, 1234 > ' ) - atomic mass= ',F7.3,' - fixed') 1235 1192 FORMAT(5X, I4, A5 ,' (',3F11.5,' )') 1236 1193 FORMAT(5X, I4, A5, ' (',3F11.5, 1237 > ' ) - atomic mass= ',F7.3,' - z fixed') 1238 1194 FORMAT(5X, I4, A5, ' (',3F11.5, 1239 > ' ) - atomic mass= ',F7.3,A) 1240 1241 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 1242 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 1243 1211 FORMAT(5X,' number of constraints = ', I6,' ( DOF = ',I6,' )' ) 1244 1219 FORMAT(/' number of active electrons: spin up=',F6.2, 1245 > ' down=',F6.2,A) 1246 1220 FORMAT(/' number of active electrons: spin up=',I6, 1247 > ' down=',I6,A) 1248 1221 FORMAT( ' number of active orbitals : spin up=',I6, 1249 > ' down=',I6,A) 1250 1251 1222 format(5x,' initial kinetic energy=',e12.5,' (ion)',2x, 1252 > e12.5,' (c.o.m.)') 1253 1223 format(5x,' after scaling= ',e12.5,' (ion)') 1254 1224 format(5x,' increased energy= ',e12.5,' (ion)') 1255 1226 format(/' final kinetic energy= ', e12.5,' (ion)',2x, 1256 > e12.5,' (c.o.m.)') 1257 1258 1230 FORMAT(/' supercell:') 1259 1231 FORMAT(5x,' omega=',F12.1) 1260 1232 FORMAT(5x,' lattice: a= ',f8.3,' b= ',f8.3,' c= ',f8.3, 1261 > /5x,' alpha=',f8.3,' beta=',f8.3,' gamma=',f8.3) 1262 1233 FORMAT(5x,' cell_name: ',A) 1263 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 1264 1242 FORMAT(5x,' a2=<',3f8.3,' >') 1265 1243 FORMAT(5x,' a3=<',3f8.3,' >') 1266 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >') 1267 1245 FORMAT(5x,' b2=<',3f8.3,' >') 1268 1246 FORMAT(5x,' b3=<',3f8.3,' >') 1269 1270 1250 FORMAT(/5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1271 & '( ',I8,' waves ',I8,' per task)') 1272 1251 FORMAT(5X,' wavefnc cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1273 & '( ',I8,' waves ',I8,' per task)') 1274 1275 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,' and',I3) 1276 1261 FORMAT(5X,' madelung=',f14.8) 1277 1270 FORMAT(/' technical parameters for minimizer:') 1278 1271 FORMAT(5x, ' translation constrained') 1279 1272 FORMAT(5x, ' rotation constrained') 1280 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1) 1281 1281 FORMAT(5X, ' maximum iterations =',I10, 1282 > ' ( ',I4,' inner ',I6,' outer )') 1283 1290 FORMAT(5X, ' tolerance=',E8.3,' (energy)',E12.3, 1284 & ' (density)') 1285 1291 FORMAT(/' Kohn-Sham scf parameters:') 1286 1292 FORMAT(5X, ' Kohn-Sham algorithm = ',A) 1287 1293 FORMAT(5X, ' SCF algorithm = ',A) 1288 1294 FORMAT(5X, ' SCF mixing parameter =',F7.4) 1289 1295 FORMAT(5X, ' Kohn-Sham iterations = ',I4, 1290 > ' (',I4,' outer)') 1291 1300 FORMAT(//) 1292 1301 FORMAT(5X, ' Kerker damping =',F7.4) 1293 1310 FORMAT(/' molecular dynamics parameters:') 1294 1320 FORMAT(5X, ' time step=',F10.2,5X,' iterations=',I10, 1295 > ' ( ',I4,' inner ',I6,' outer )') 1296 1321 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1, 1297 > /5X,' iterations=',I10, 1298 > ' ( ',I4,' inner ',I6,' outer )') 1299 1330 FORMAT(5X, ' integration algorithm= ',A) 1300 1340 FORMAT(/5X, ' cooling/heatting rate= ',e12.5,' (ion)') 1301 1350 FORMAT(I8,2E19.10,E14.5,4F14.2,2E19.10,2F8.4) 1302 1360 format(100e19.10) 1303 1304 1393 format(/' Pressure Output Generated ') 1305 1394 format(/' Constant Energy Simulation ') 1306 1395 format(/' Nose-Hoover Simulation - Thermostat Parameters:') 1307 1396 format(5x, 'SA decay rates =',e10.3,' (ion)') 1308 1397 format(5x, 'link = ',I3, 1309 > ' Te =',f8.2,' Qe =',e10.3,' 2*pi/we=',e10.3,' Ee0=',e10.3) 1310 1398 format(5x, 'link = ',I3, 1311 > ' Tr =',f8.2,' Qr =',e10.3,' 2*pi/wr=',e10.3,' Er0=',e10.3) 1312 1313 1399 format(//' Mulliken Analysis Output Generated ') 1314 1400 FORMAT(I3,3E18.8/3X,3E18.8) 1315 1410 FORMAT(10X,'============= summary of results =================') 1316 1420 FORMAT(/' final position of ions:') 1317 1421 FORMAT(/' final velocity of ions:') 1318 1471 FORMAT(/' Kinetic energy (elc) :',E19.10,' (',E15.5,'/elc)') 1319 1472 FORMAT( ' Kinetic energy (ion) :',E19.10,' (',E15.5,'/ion)') 1320 1473 FORMAT( ' thermostat energy (elc) :',E19.10,' (',E15.5,'/elc)') 1321 1474 FORMAT( ' thermostat energy (ion) :',E19.10,' (',E15.5,'/ion)') 1322 1480 FORMAT(' Temperature : ',F10.1,' K (ion)') 1323 1490 FORMAT(' : ',F10.1,' K (c.o.m.)') 1324 1491 FORMAT(' Temperature : ',F10.1,' K (elc)') 1325 1492 FORMAT(/' Eaverage : ',E19.10) 1326 1493 FORMAT( ' Evariance : ',E19.10) 1327 1494 FORMAT( ' Cv - f*kb/(2*nion) : ',E19.10) 1328 1600 FORMAT(/' Total PSPW energy :',E19.10) 1329 1801 FORMAT(//'== Timing ==') 1330 1802 FORMAT(/'cputime in seconds') 1331 1803 FORMAT( ' prologue : ',E14.6) 1332 1804 FORMAT( ' main loop : ',E14.6) 1333 1805 FORMAT( ' epilogue : ',E14.6) 1334 1806 FORMAT( ' total : ',E14.6) 1335 1807 FORMAT( ' cputime/step: ',E14.6, 1336 > ' (',I8,' evalulations,', I8,' linesearches)') 1337 1808 FORMAT(A,E14.6,E14.6) 1338 1809 FORMAT(//A,2A14) 1339 1340 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 1341 1342 9000 if (taskid.eq.MASTER) write(6,9010) ierr 1343 call Parallel_Finalize() 1344 1345 cgmdv5 = value 1346 return 1347 END 1348