1* 2* $Id$ 3* 4 5*********************************************************************** 6* band_sd * 7* * 8* This is a developing steepest descent code for BAND * 9* * 10* * 11*********************************************************************** 12 13 logical function band_sd(rtdb) 14 implicit none 15 integer rtdb 16 17#include "global.fh" 18#include "bafdecls.fh" 19#include "inp.fh" 20#include "btdb.fh" 21#include "stdio.fh" 22#include "util.fh" 23#include "errquit.fh" 24 25 26* **** parallel variables **** 27 integer taskid,np,np_i,np_j,np_k 28 integer MASTER 29 parameter(MASTER=0) 30 31* **** timing variables **** 32 real*8 cpu1,cpu2,cpu3,cpu4 33 real*8 t1,t2,t3,t4,av 34 35* **** lattice variables **** 36 integer ngrid(3),nwave,nfft3d 37 integer npack1,npack0 38 39* **** electronic variables **** 40 logical spin_orbit,newpsi 41 real*8 icharge 42 integer ispin,ispinq,nbrillioun,nbrillq 43 integer ne(2),n1(2),n2(2),nemax,neq(2),nemaxq,neall 44 real*8 en(2) 45 real*8 dipole(3) 46 47 integer psi1_tag,psi2_tag,psir_tag,Hpsi_tag,next 48 integer psi1_shift,psi2_shift 49 integer dn(2) 50 51 52* ***** energy variables **** 53 real*8 E(50) 54 55 integer eig_tag,hml_tag,lmd_tag,svec_tag 56 integer eig_shift,hml_shift,lmd_shift,svec_shift 57 58* **** psi smearing block **** 59 logical fractional 60 integer smearoccupation,smeartype 61 real*8 smearfermi(2),smearcorrection,smearkT 62 63 64* **** error variables **** 65 integer ierr 66 67* **** local variables **** 68 complex*16 zero,one 69 parameter (zero=(0.0d0,0.0d0), one=(1.0d0,0.0d0)) 70 71 logical lprint,mprint,hprint 72 integer ms,mapping,mapping1d 73 real*8 deltae,deltac,deltar 74 real*8 gx,gy,gz,cx,cy,cz,sum1,sum2 75 real*8 EV,pi,f0,f1,f2,f3,f4,f5,f6 76 integer i,j,k,ia,n,nn,nb 77 integer ii,jj,indx,indx1,if1,if2 78 integer icount,it_in,it_out 79 real*8 w,sumall,virial,a,b,c,alpha,beta,gamma 80 integer nfft3 81 parameter (nfft3=32) 82 character*255 full_filename 83 84 logical value,psi_nogrid 85 integer hversion,hnfft(3),hispin,hne(2) 86 real*8 hunita(3,3) 87 integer ind,vers 88 89 character*255 filename 90 character*50 control_input_psi 91 external control_input_psi 92 logical c_wvfnc_expander 93 external c_wvfnc_expander 94 95 96* **** external functions **** 97 real*8 cpsp_zv,cpsp_rc,ewald_rcut,ion_amass,ion_TotalCharge 98 real*8 ewald_mandelung 99 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 100 real*8 lattice_unitg,brillioun_weight 101 integer ewald_ncut,ewald_nshl3d 102 integer cpsp_nprj,cpsp_lmax,cpsp_locp,cpsp_psp_type 103 character*4 ion_aname,ion_atom 104 external cpsp_zv,cpsp_rc,ewald_rcut,ion_amass,ion_TotalCharge 105 external ewald_mandelung 106 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 107 external lattice_unitg,brillioun_weight 108 external ewald_ncut,ewald_nshl3d 109 external cpsp_nprj,cpsp_lmax,cpsp_locp,cpsp_psp_type 110 external ion_aname,ion_atom 111 112 logical control_fractional 113 external control_fractional 114 real*8 control_fractional_temperature 115 external control_fractional_temperature 116 real*8 control_fractional_kT,control_fractional_smeartype 117 external control_fractional_kT,control_fractional_smeartype 118 real*8 control_tole,control_tolc,control_tolr,ion_rion 119 external control_tole,control_tolc,control_tolr,ion_rion 120 real*8 control_time_step,control_fake_mass 121 external control_time_step,control_fake_mass 122 logical control_read,control_move,ion_init,ion_q_FixIon 123 external control_read,control_move,ion_init,ion_q_FixIon 124 logical ion_q_xyzFixIon,band_HFX,band_HFX_relaxed 125 external ion_q_xyzFixIon,band_HFX,band_HFX_relaxed 126 character*14 ion_q_xyzFixIon_label 127 external ion_q_xyzFixIon_label 128 129 logical brillioun_print 130 external brillioun_print 131 real*8 brillioun_k_brdcst,brillioun_ks_brdcst 132 external brillioun_k_brdcst,brillioun_ks_brdcst 133 real*8 brillioun_weight_brdcst 134 external brillioun_weight_brdcst 135 real*8 cpsi_eig_brdcst_tag,cpsi_occ_brdcst_tag 136 external cpsi_eig_brdcst_tag,cpsi_occ_brdcst_tag 137 real*8 cpsi_sv_brdcst_tag 138 external cpsi_sv_brdcst_tag 139 140 integer Cram_nwave_brdcst,Cram_nwave_all_brdcst 141 external Cram_nwave_brdcst,Cram_nwave_all_brdcst 142 143 integer psi_get_version,brillioun_nbrillioun 144 integer Pneb_ispinq,Pneb_nbrillq,Pneb_w_size 145 integer cpsi_data_alloc,cpsi_data_get_next,cpsi_data_get_chnk 146 integer control_it_in,control_it_out,control_gga,control_version 147 integer control_ngrid,pack_nwave 148 integer ion_nion,ion_natm,ion_katm,ion_nkatm 149 external psi_get_version,brillioun_nbrillioun 150 external Pneb_ispinq,Pneb_nbrillq,Pneb_w_size 151 external cpsi_data_alloc,cpsi_data_get_next,cpsi_data_get_chnk 152 external control_it_in,control_it_out,control_gga,control_version 153 external control_ngrid,pack_nwave 154 external ion_nion,ion_natm,ion_katm,ion_nkatm 155 156 character*12 control_boundry 157 external control_boundry 158 159 logical pspw_reformat_c_wvfnc 160 logical pspw_HFX,pspw_HFX_relaxed 161 logical cpsp_semicore,control_Mulliken 162 real*8 cpsp_rcore,cpsp_ncore 163 external pspw_reformat_c_wvfnc 164 external pspw_HFX,pspw_HFX_relaxed 165 external cpsp_semicore,control_Mulliken 166 external cpsp_rcore,cpsp_ncore 167 logical control_check_charge_multiplicity 168 external control_check_charge_multiplicity 169 real*8 nwpw_timing 170 external nwpw_timing 171 integer control_np_dimensions 172 integer control_mapping,control_mapping1d 173 external control_np_dimensions 174 external control_mapping,control_mapping1d 175 176 logical control_print 177 logical control_translation,control_rotation,control_balance 178 external control_print 179 external control_translation,control_rotation,control_balance 180 181 logical Dneall_m_allocate,Dneall_m_free,ion_disp_on 182 external Dneall_m_allocate,Dneall_m_free,ion_disp_on 183 184 complex*16 Pneb_w_value 185 external Pneb_w_value 186 character*9 ion_amm 187 external ion_amm 188 189 character*255 cpsp_comment,comment 190 external cpsp_comment 191 integer ion_nconstraints,ion_ndof 192 external ion_nconstraints,ion_ndof 193 logical ion_makehmass2 194 external ion_makehmass2 195 196 197 198* |************| 199*****************************| PROLOGUE |**************************** 200* |************| 201 202 value = .true. 203 pi = 4.0d0*datan(1.0d0) 204 205 call nwpw_timing_init() 206 call dcopy(30,0.0d0,0,E,1) 207 208 209* **** get parallel variables **** 210 call Parallel_Init() 211 call Parallel_np(np) 212 call Parallel_taskid(taskid) 213 214 if (.not.control_read(13,rtdb)) 215 > call errquit('band_sd:error reading control',0,DISK_ERR) 216 217 218 lprint = ((taskid.eq.MASTER).and.(control_print(print_low))) 219 mprint = ((taskid.eq.MASTER).and.(control_print(print_medium))) 220 hprint = ((taskid.eq.MASTER).and.(control_print(print_high))) 221 222 if (taskid.eq.MASTER) call current_second(cpu1) 223 224* ***** print out header **** 225 if (mprint) then 226 write(luout,1000) 227 write(luout,1010) 228 write(luout,1020) 229 write(luout,1010) 230 write(luout,1030) 231 write(luout,1010) 232 write(luout,1035) 233 write(luout,1010) 234 write(luout,1040) 235 write(luout,1010) 236 write(luout,1041) 237 write(luout,1042) 238 write(luout,1043) 239 write(luout,1010) 240 write(luout,1000) 241 call nwpw_message(1) 242 write(luout,1110) 243 end if 244 245 call Parallel3d_Init(control_np_dimensions(2), 246 > control_np_dimensions(3)) 247 call Parallel3d_np_i(np_i) 248 call Parallel3d_np_j(np_j) 249 call Parallel3d_np_k(np_k) 250 251 ngrid(1) = control_ngrid(1) 252 ngrid(2) = control_ngrid(2) 253 ngrid(3) = control_ngrid(3) 254 nwave = 0 255 mapping = control_mapping() 256 257* **** initialize D3dB data structure **** 258 call C3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 259 call C3dB_nfft3d(1,nfft3d) 260 261* **** initialize psi_data **** 262 call cpsi_data_init(20) 263 264* **** read ions **** 265 value = ion_init(rtdb) 266 call center_geom(cx,cy,cz) 267 call center_mass(gx,gy,gz) 268 269* **** initialize lattice and packing data structure **** 270 call lattice_init() 271 call c_G_init() 272 call brillioun_init() 273 call Cram_Init() 274 call C3dB_pfft_init() 275 276* ***** Initialize double D3dB data structure **** 277 if ((control_gga().ge.10).and.(control_gga().le.200)) then 278 call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 279 call G_init() 280 call mask_init() 281 end if 282 283c* **** read ions **** 284c value = ion_init(rtdb) 285c call center_geom(cx,cy,cz) 286c call center_mass(gx,gy,gz) 287 288* **** allocate psp data structure and read in psedupotentials into it **** 289 call cpsp_init() 290 call cpsp_readall() 291 if (cpsp_semicore(0)) call c_semicore_check() 292 293 294* **** initialize G,mask,ke,and coulomb data structures **** 295 call cstrfac_init() 296 call cke_init() 297 call c_coulomb_init() 298 call ewald_init() 299 300 301* **** generate initial wavefunction if it does not exist **** 302 if (.not.control_check_charge_multiplicity()) then 303 call cpsi_new() 304 newpsi = .true. 305 else 306 newpsi = .false. 307 308* **** convert from pspw format to band format **** 309 vers = psi_get_version() 310 if ((vers.eq.3).or.(vers.eq.4)) then 311 newpsi = .true. 312 value = btdb_parallel(.false.) 313 if (taskid.eq.MASTER) then 314 value= pspw_reformat_c_wvfnc(1) 315 end if 316 value = btdb_parallel(.true.) 317 end if 318 end if 319 320 call psi_get_ne(ispin,ne) 321 if (ispin.eq.3) then 322 spin_orbit = .true. 323 ispin=2 324 else 325 spin_orbit = .false. 326 end if 327 nbrillioun = brillioun_nbrillioun() 328 call Pneb_init(ispin,ne,nbrillioun,spin_orbit) 329 call Pneb_neq(neq) 330 331 332* ***** allocate psi2,and psi1 wavefunctions **** 333 call psi_get_ne_occupation(ispin,ne,smearoccupation) 334 if (smearoccupation.gt.0) then 335 fractional = .true. 336 else 337 fractional = .false. 338 end if 339 mapping1d = control_mapping1d() 340 ispinq = Pneb_ispinq() 341 nbrillq = Pneb_nbrillq() 342 343 call Cram_npack(0,npack0) 344 call Cram_max_npack(npack1) 345 nemaxq = neq(1)+neq(2) 346 neall = ne(1) +ne(2) 347 348 psi2_tag = cpsi_data_alloc(nbrillq,nemaxq,2*npack1) 349 psi1_tag = cpsi_data_alloc(nbrillq,nemaxq,2*npack1) 350 351* *** fractional orbitals *** 352 if (smearoccupation.gt.0) then 353 call cpsi_data_set_next(psi2_tag, 354 > cpsi_data_alloc(nbrillq,nemaxq,1)) 355 call cpsi_data_set_next(psi1_tag, 356 > cpsi_data_alloc(nbrillq,nemaxq,1)) 357 smeartype = control_fractional_smeartype() 358 smearkT = control_fractional_kT() 359 end if 360 361c **** allocate other variables **** 362 Hpsi_tag = cpsi_data_alloc(nbrillq,nemaxq,2*npack1) 363 psir_tag = cpsi_data_alloc(nbrillq,nemaxq,2*nfft3d) 364 hml_tag = cpsi_data_alloc(nbrillq,1,2*Pneb_w_size(0,1)) 365 eig_tag = cpsi_data_alloc(nbrillq,ne(1)+ne(2),1) 366 svec_tag = cpsi_data_alloc(nbrillq,neq(1),3) 367 value = BA_alloc_get(mt_dbl,2*nfft3d,'dn',dn(2),dn(1)) 368 if (.not. value) 369 > call errquit('band_sd:out of heap memory',0,MA_ERR) 370 371 372* ***** read initial wavefunctions into psi2 **** 373 if (.not.btdb_get(rtdb,'nwpw:psi_nogrid', 374 > mt_log,1,psi_nogrid)) 375 > psi_nogrid = .true. 376 377 if (psi_nogrid) then 378 379 call psi_get_header(hversion,hnfft,hunita,hispin,hne) 380 381 if ( (hnfft(1).ne.control_ngrid(1)) .or. 382 > (hnfft(2).ne.control_ngrid(2)) .or. 383 > (hnfft(3).ne.control_ngrid(3)) ) then 384 385 hnfft(1) = control_ngrid(1) 386 hnfft(2) = control_ngrid(2) 387 hnfft(3) = control_ngrid(3) 388 call Parallel_taskid(taskid) 389 390 call ga_sync() 391 value = btdb_parallel(.false.) 392 call ga_sync() 393 if (hprint) then 394 395 filename = control_input_psi() 396 397 ind = index(filename,' ') - 1 398 if (.not. btdb_cput(rtdb,'c_xpndr:old_wavefunction_filename', 399 > 1,filename(1:ind))) 400 > call errquit( 401 > 'c_wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 402 403 if (.not. btdb_cput(rtdb,'c_xpndr:new_wavefunction_filename', 404 > 1,filename(1:ind))) 405 > call errquit( 406 > 'c_wvfnc_expander_input: btdb_cput failed', 0, RTDB_ERR) 407 408 if (.not. btdb_put(rtdb,'c_xpndr:ngrid',mt_int,3,hnfft)) 409 > call errquit( 410 > 'c_wvfnc_expander_input: btdb_put failed', 0, RTDB_ERR) 411 412 write(*,*) 413 write(*,*) "Grid is being converted:" 414 write(*,*) "------------------------" 415 write(*,*) 416 write(*,*) "To turn off automatic grid conversion:" 417 write(*,*) 418 write(*,*) "set nwpw:psi_nogrid .false." 419 write(*,*) 420 value = c_wvfnc_expander(rtdb) 421 422 end if 423 call ga_sync() 424 value = btdb_parallel(.true.) 425 426 end if 427 end if 428 429* ***** read psi2 wavefunctions **** 430 call cpsi_read(spin_orbit,ispin,ne,nbrillioun,psi2_tag) 431 432* **** Ortho Check **** 433 call Pneb_orthoCheckMake_tag(.true.,0,0,npack1,psi2_tag) 434 435* **** initialize SIC and HFX **** 436c call pspw_init_SIC(rtdb,ne) 437 call band_init_HFX(rtdb,nbrillioun,ispin,ne) 438 439* **** initialize QM/MM **** 440c call pspw_init_APC(rtdb) 441c call pspw_qmmm_init(rtdb) 442 443 444* **** initialize FixIon constraint **** 445 call ion_init_FixIon(rtdb) 446 447 448* |**************************| 449****************** summary of input data ********************** 450* |**************************| 451 452 453* **** determine en **** 454 if (.not.spin_orbit) then 455 next = cpsi_data_get_next(psi2_tag) 456 icharge = 0.0d0 457 en(1) = 0.0d0 458 en(2) = 0.0d0 459 b = dble(3-ispin) 460 do nb=1,nbrillq 461 w = brillioun_weight(nb) 462 do ms=1,ispin 463 do i=1,ne(ms) 464 if (next.lt.0) then 465 a = 1.0d0 466 else 467 indx = cpsi_data_get_chnk(next,nb)+i-1 468 if (.not.spin_orbit) indx = indx + (ms-1)*ne(1) 469 a = dbl_mb(indx) 470 end if 471 icharge = icharge - b*a*w 472 en(ms) = en(ms) + a*w 473 end do 474 end do 475 end do 476 call K1dB_Vector_SumAll(2,en) 477 call K1dB_SumAll(icharge) 478 else 479 icharge = -ne(1) 480 en(1) = ne(1) 481 en(ispin) = ne(ispin) 482 end if 483 484 485 if (mprint) then 486 write(luout,1111) np 487 write(luout,1117) np_i,np_j,np_k 488 if (mapping.eq.1) write(luout,1112) 489 if (mapping.eq.2) write(luout,1113) 490 if (mapping.eq.3) write(luout,1118) 491 492 write(luout,1115) 493 IF(control_move()) THEN 494 write(luout,1120) 'yes' 495 ELSE 496 write(luout,1120) 'no' 497 ENDIF 498 write(luout,1121) control_boundry(),control_version() 499 if (.not.spin_orbit) then 500 if (ispin.eq.1) write(luout,1130) 'restricted' 501 if (ispin.eq.2) write(luout,1130) 'unrestricted' 502 else 503 write(luout,1130) 'spin orbit' 504 end if 505 506 call v_bwexc_print(luout,control_gga()) 507 508c if (fractional) write(6,1132) 509c call pspw_print_SIC(6) 510c call pspw_print_HFX(6) 511 if (ion_makehmass2()) write(luout,1135) 512 write(luout,1140) 513 do ia = 1,ion_nkatm() 514 write(luout,1150) ia,ion_atom(ia), 515 > cpsp_zv(ia),cpsp_lmax(ia) 516 517 comment = cpsp_comment(ia) 518 i = inp_strlen(comment) 519 write(luout,1157) comment(1:i) 520 write(luout,1158) cpsp_psp_type(ia) 521 write(luout,1152) cpsp_lmax(ia) 522 write(luout,1153) cpsp_locp(ia) 523 write(luout,1154) cpsp_nprj(ia) 524 if (cpsp_semicore(ia)) 525 > write(luout,1155) cpsp_rcore(ia),cpsp_ncore(ia) 526 write(luout,1151) (cpsp_rc(i,ia),i=0,cpsp_lmax(ia)) 527 end do 528 529 530 icharge = icharge + ion_TotalCharge() 531 write(6,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),(ion_rion(K,I),K=1,3), 539 > ion_amass(I)/1822.89d0,ion_amm(i) 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),(ion_rion(K,I),K=1,3), 545 > ion_amass(I)/1822.89d0,ion_amm(i) 546 end if 547 end do 548 write(luout,1200) cx,cy,cz 549 write(luout,1210) gx,gy,gz 550 write(luout,1211) ion_nconstraints(),ion_ndof() 551 552 write(luout,1220) en(1),en(ispin),' (Fourier space)' 553 write(luout,1221) ne(1),neq(1), 554 > ne(ispin),neq(ispin),' (Fourier space)' 555 556 write(luout,1230) 557 write(luout,1241) lattice_unita(1,1), 558 > lattice_unita(2,1), 559 > lattice_unita(3,1) 560 write(luout,1242) lattice_unita(1,2), 561 > lattice_unita(2,2), 562 > lattice_unita(3,2) 563 write(luout,1243) lattice_unita(1,3), 564 > lattice_unita(2,3), 565 > lattice_unita(3,3) 566 write(luout,1244) lattice_unitg(1,1), 567 > lattice_unitg(2,1), 568 > lattice_unitg(3,1) 569 write(luout,1245) lattice_unitg(1,2), 570 > lattice_unitg(2,2), 571 > lattice_unitg(3,2) 572 write(luout,1246) lattice_unitg(1,3), 573 > lattice_unitg(2,3), 574 > lattice_unitg(3,3) 575 write(luout,1231) lattice_omega() 576 call lattice_abc_abg(a,b,c,alpha,beta,gamma) 577 write(luout,1232) a,b,c,alpha,beta,gamma 578 write(luout,1260) ewald_rcut(),ewald_ncut() 579 write(luout,1261) ewald_mandelung() 580 581 write(luout,1255) 582 write(luout,1256) brillioun_nbrillioun() 583 end if 584 585 586c write(6,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 587c > pack_nwave_all(0),pack_nwave(0) 588c write(6,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 589c > pack_nwave_all(1),pack_nwave(1) 590 591c **** print brillioun zone - extra logic for distributed kpoints **** 592 if (brillioun_print()) then 593 do i=1,brillioun_nbrillioun() 594 f0 = brillioun_weight_brdcst(i) 595 f1 = brillioun_ks_brdcst(1,i) 596 f2 = brillioun_ks_brdcst(2,i) 597 f3 = brillioun_ks_brdcst(3,i) 598 f4 = brillioun_k_brdcst(1,i) 599 f5 = brillioun_k_brdcst(2,i) 600 f6 = brillioun_k_brdcst(3,i) 601 if (mprint) write(luout,1257) f0,f1,f2,f3,f4,f5,f6 602 end do 603 else 604 if (mprint) write(luout,1258) 605 end if 606 607 if1 = Cram_nwave_all_brdcst(0) 608 if2 = Cram_nwave_brdcst(0) 609 if (mprint) then 610 write(luout,1249) 611 write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 612 > if1,if2 613 end if 614 615 if (brillioun_print()) then 616 do i=1,brillioun_nbrillioun() 617 if1 = Cram_nwave_all_brdcst(i) 618 if2 = Cram_nwave_brdcst(i) 619 if (mprint) then 620 write(luout,1251) i,lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 621 > if1,if2 622 end if 623 end do 624 else 625 if (mprint) write(luout,1252) 626 end if 627 628 629 if (mprint) then 630 631 write(luout,1270) 632 if (.not.control_translation()) write(luout,1271) 633 if (.not.control_rotation()) write(luout,1272) 634 write(luout,1280) control_time_step(),control_fake_mass() 635 write(luout,1290) control_tole(),control_tolc(),control_tolr() 636 write(luout,1281) control_it_in()*control_it_out(), 637 > control_it_in(),control_it_out() 638 639 if (control_fractional()) then 640 write(6,1297) 641 if (control_fractional_smeartype().eq.0) 642 > write(6,1298) "step function" 643 if (control_fractional_smeartype().eq.1) 644 > write(6,1298) "Fermi-Dirac" 645 if (control_fractional_smeartype().eq.2) 646 > write(6,1298) "Gaussian" 647 write(6,1299) control_fractional_kT(), 648 > control_fractional_temperature() 649 end if 650 write(luout,1300) 651 write(luout,1305) 652 call util_flush(luout) 653 end if 654c 655c* |***************************| 656c****************** start iterations ********************** 657c* |***************************| 658c 659 if (taskid.eq.MASTER) call current_second(cpu2) 660 if (taskid.eq.MASTER) CALL nwpw_MESSAGE(2) 661 it_in = control_it_in() 662 it_out = control_it_out() 663 icount = 0 664 1 continue 665 icount = icount + 1 666 call band_inner_loop(ispin,ispinq,ne,neq,nbrillioun,nbrillq, 667 > nfft3d, 668 > psi1_tag,psi2_tag,dbl_mb(dn(1)), 669 > it_in,E,deltae,deltac,deltar, 670 > hml_tag, 671 > psir_tag,Hpsi_tag) 672 673 if (mprint) then 674 write(luout,1310) icount*it_in,E(1),deltae,deltac,deltar 675 call util_flush(luout) 676 end if 677 if ((deltae.gt.0.0d0).and.(icount.gt.1)) then 678 if (mprint) 679 > write(luout,*) 680 > ' *** Energy going up. iteration terminated.' 681 go to 2 682 end if 683 deltae = dabs(deltae) 684 if ((deltae.lt.control_tole()).and. 685 > (deltac.lt.control_tolc()).and. 686 > (deltar.lt.control_tolr())) then 687 if (mprint) 688 > write(luout,*) 689 > ' *** tolerance ok. iteration terminated.' 690 go to 2 691 end if 692 if (icount.lt.it_out) go to 1 693 if (mprint) 694 > write(luout,*) 695 > '*** arived at the Maximum iiteration. terminated.' 696 697*:::::::::::::::::::: end of iteration loop ::::::::::::::::::::::::: 698 699 2 continue 700 if (taskid.eq.MASTER) CALL nwpw_MESSAGE(3) 701 if (taskid.eq.MASTER) call current_second(cpu3) 702 703 704 705* |****************************************| 706*********** produce CHECK file and diagonalize hml ***************** 707* |****************************************| 708 709* **** produce CHECK FILE **** 710 if (taskid.eq.MASTER) then 711 call util_file_name('CHECK',.true., 712 > .false., 713 > full_filename) 714 open(unit=17,file=full_filename,form='formatted') 715 end if 716 717* **** check total number of electrons **** 718 do ms =1,ispin 719 call C3dB_r_dsum(1,dbl_mb(dn(1)+(ms-1)*nfft3d),sumall) 720 en(ms) = sumall*lattice_omega() 721 > /dble(ngrid(1)*ngrid(2)*ngrid(3)) 722 end do 723 if (taskid.eq.MASTER) then 724 write(17,1320) (en(ms),ms=1,ispin) 725 end if 726 727* **** comparison between hamiltonian an lambda matrix **** 728* **** not done - because this can generate a very large data file **** 729 730* **** check orthonormality **** 731 if (taskid.eq.MASTER) then 732 write(17,1350) 733 end if 734 735c lmd_tag = cpsi_data_alloc(1,1,2*Pneb_w_size(0,1)) 736c lmd_shift = cpsi_data_get_chnk(lmd_tag,1) 737c do nb=1,nbrillq 738c psi1_shift = cpsi_data_get_chnk(psi1_tag,nb) 739c call Pneb_ffw_Multiply(nb,0,dbl_mb(psi1_shift), 740c > dbl_mb(psi1_shift), 741c > dbl_mb(lmd_shift)) 742c do ms=1,ispin 743c do j=1,ne(ms) 744c do i=j,ne(ms) 745c if (taskid.eq.MASTER) 746c > write(17,1360) ms,i,j, 747c > Pneb_w_value(nb,ms,i,j,dbl_mb(lmd_shift)) 748c end do 749c end do 750c end do 751c end do 752c call cpsi_data_dealloc(lmd_tag) 753 754* **** close check file **** 755 if (taskid.eq.MASTER) then 756 close(17) 757 end if 758 759 760 761* ***** diagonalize the hamiltonian matrix but don't rotate **** 762 call cpsi_data_update(hml_tag) 763 call cpsi_data_update(eig_tag) 764 do nb=1,nbrillq 765 hml_shift = cpsi_data_get_chnk(hml_tag,nb) 766 eig_shift = cpsi_data_get_chnk(eig_tag,nb) 767 call Pneb_w_diag(0,nb,dbl_mb(eig_shift),dbl_mb(hml_shift)) 768 end do 769 call cpsi_data_noupdate(hml_tag) 770 call cpsi_data_noupdate(eig_tag) 771 772 773* ***** diagonalize and rotate the hamiltonian matrix **** 774 call cpsi_data_update(psi2_tag) 775 do nb=1,nbrillq 776 psi1_shift = cpsi_data_get_chnk(psi1_tag,nb) 777 psi2_shift = cpsi_data_get_chnk(psi2_tag,nb) 778 hml_shift = cpsi_data_get_chnk(hml_tag,nb) 779 call Pneb_fwf_Multiply(0,nb, 780 > one, 781 > dbl_mb(psi1_shift),npack1, 782 > dbl_mb(hml_shift), 783 > zero, 784 > dbl_mb(psi2_shift)) 785 end do 786 call cpsi_data_noupdate(psi2_tag) 787 788 789 790* |***************************| 791****************** report summary of results ********************** 792* |***************************| 793 call center_geom(cx,cy,cz) 794 call center_mass(gx,gy,gz) 795 796 if (mprint) then 797 write(luout,1300) 798 write(luout,1410) 799 write(luout,1420) 800 do I=1,ion_nion() 801 if (ion_q_FixIon(I)) then 802 write(luout,1191) I,ion_aname(I),(ion_rion(k,i),K=1,3), 803 > ion_amass(I)/1822.89d0,ion_amm(i) 804 else if (ion_q_xyzFixIon(I)) then 805 write(6,1194) I,ion_aname(I),(ion_rion(k,i),K=1,3), 806 > ion_amass(I)/1822.89d0,ion_q_xyzFixIon_label(I) 807 else 808 write(luout,1190) I,ion_aname(I),(ion_rion(k,i),K=1,3), 809 > ion_amass(I)/1822.89d0,ion_amm(i) 810 end if 811 end do 812 write(luout,1200) cx,cy,cz 813 write(luout,1210) gx,gy,gz 814 write(luout,1211) ion_nconstraints(),ion_ndof() 815 816 817 write(luout,*) 818 write(luout,1320) en(1),en(ispin),' (real space)' 819 820 write(luout,1430) E(1),E(1)/ion_nion() 821 write(luout,1440) E(2),E(2)/neall 822 write(luout,1450) E(3),E(3)/neall 823 write(luout,1460) E(4),E(4)/neall 824 if (band_HFX()) then 825 write(luout,1457) E(20),E(20)/n2(ispin) 826 end if 827 828 write(luout,1470) E(5),E(5)/ion_nion() 829 write(luout,1480) E(6),E(6)/neall 830 write(luout,1490) E(7),E(7)/neall 831 write(luout,1495) E(8),E(8)/neall 832 write(luout,1496) E(9),E(9)/neall 833 write(luout,1497) E(10),E(10)/neall 834 if (band_HFX().and.band_HFX_relaxed()) then 835 write(luout,1502) E(21),E(21)/n2(ispin) 836 end if 837 virial = (E(10)+E(9)+E(8)+E(7))/E(6) 838 write(luout,1498) virial 839 840 if (ion_disp_on()) then 841 write(luout,1720) E(33) 842 end if 843 844 end if 845 846 NN=ne(1)-ne(2) 847 EV=27.2116d0 848 if (mprint) then 849 if (control_fractional()) then 850 if (ispin.eq.1) then 851 write(luout,1507) smearfermi(1),smearfermi(1)*EV 852 else 853 write(luout,1507) smearfermi(1),smearfermi(1)*EV, 854 > smearfermi(2),smearfermi(2)*EV 855 end if 856 end if 857 end if 858 859* *** generate spinorbit vector *** 860 if (spin_orbit) call Pneb_f_SOSpins_tag(psi2_tag,svec_tag) 861 862* *** printout the eigenvalue spetra **** 863 if (brillioun_print()) then 864 do nb=1,brillioun_nbrillioun() 865 f0 = brillioun_weight_brdcst(nb) 866 f1 = brillioun_ks_brdcst(1,nb) 867 f2 = brillioun_ks_brdcst(2,nb) 868 f3 = brillioun_ks_brdcst(3,nb) 869 f4 = brillioun_k_brdcst(1,nb) 870 f5 = brillioun_k_brdcst(2,nb) 871 f6 = brillioun_k_brdcst(3,nb) 872 if (mprint) then 873 write(luout,1508) nb,f0,f1,f2,f3,f4,f5,f6 874 write(luout,1500) 875 end if 876 next = cpsi_data_get_next(psi1_tag) 877 if (spin_orbit) then 878c if (mprint) write(luout,1511) 879 do i=0,ne(1)-1 880 f1=cpsi_eig_brdcst_tag(ne,spin_orbit,eig_tag, nb,1,ne(1)-i) 881 f2=cpsi_sv_brdcst_tag( spin_orbit,svec_tag,nb,ne(1)-i,1) 882 f3=cpsi_sv_brdcst_tag( spin_orbit,svec_tag,nb,ne(1)-i,2) 883 f4=cpsi_sv_brdcst_tag( spin_orbit,svec_tag,nb,ne(1)-i,3) 884 f0 = dsqrt(f2*f2 + f3*f3 + f4*f4) 885 f5=cpsi_occ_brdcst_tag(ne,spin_orbit,next, nb,1,ne(1)-i) 886 if (mprint) write(luout,1512) f1,f1*EV,f0,f2,f3,f4,f5 887 end do 888 else 889 do i=0,NN-1 890 f1=cpsi_eig_brdcst_tag(ne,spin_orbit,eig_tag,nb,1,ne(1)-i) 891 f2=cpsi_occ_brdcst_tag(ne,spin_orbit,next, nb,1,ne(1)-i) 892 if (mprint) write(luout,1510) f1,f1*EV,f2 893 end do 894 do i=0,ne(2)-1 895 f1=cpsi_eig_brdcst_tag(ne,spin_orbit,eig_tag,nb,1,ne(1)-i-NN) 896 f2=cpsi_occ_brdcst_tag(ne,spin_orbit,next, nb,1,ne(1)-i-NN) 897 f3=cpsi_eig_brdcst_tag(ne,spin_orbit,eig_tag,nb,2,ne(2)-i) 898 f4=cpsi_occ_brdcst_tag(ne,spin_orbit,next, nb,2,ne(2)-i) 899 if (mprint) write(luout,1510) f1,f1*EV,f2,f3,f3*EV,f4 900 end do 901 end if 902 end do 903 endif 904 905 906 907* ***** extra energy output for QA test **** 908 if (mprint) then 909 write(luout,1600) E(1) 910 end if 911 912* |***************************| 913****************** Prologue ********************** 914* |***************************| 915c 916c* **** calculate spin contamination **** 917c call Calculate_psi_spin2(ispin,ne,npack1,dcpl_mb(psi2(1)), 918c > fractional,dbl_mb(occ2(1)),w) 919c 920c* **** calculate the Dipole *** 921c call Calculate_Dipole(ispin,ne,n2ft3d,dbl_mb(dn(1)),dipole) 922c 923c* **** perfom Lubin and Mulliken analysis *** 924c if (control_Mulliken()) then 925c 926c* **** Lubin Water Analysis *** 927c call pspw_Lubin_water_analysis(rtdb,ispin,ne,n2ft3d, 928c > dbl_mb(dn(1))) 929c 930c* **** Analysis *** 931c call pspw_analysis(0,rtdb,ispin,ne,dcpl_mb(psi2(1)), 932c > dbl_mb(eig(1))) 933c 934c* **** generate APC ***** 935c call pspw_dngen_APC(ispin,ne,dbl_mb(dn(1))) 936c call pspw_print_APC(6) 937c 938c end if 939 940* ***** write psi2 wavefunctions **** 941 call cpsi_write(spin_orbit,ispin,ne,nbrillioun,psi2_tag) 942 943* **** write geometry to rtdb **** 944 call ion_write(rtdb) 945 946 947* **** deallocate heap memory **** 948 call ewald_end() 949 call cstrfac_end() 950 call c_coulomb_end() 951 call cke_end() 952 call cpsp_end() 953 call Cram_end() 954 call c_G_end() 955 call brillioun_end() 956 957 call ion_end() 958 call ion_end_FixIon() 959c call pspw_end_SIC() 960 call band_end_HFX() 961c call pspw_end_APC() 962c call pspw_qmmm_end() 963c 964 call cpsi_data_dealloc(Hpsi_tag) 965 call cpsi_data_dealloc(psir_tag) 966 next = cpsi_data_get_next(psi1_tag) 967 if (next.ge.0) call cpsi_data_dealloc(next) 968 call cpsi_data_dealloc(psi1_tag) 969 970 next = cpsi_data_get_next(psi2_tag) 971 if (next.ge.0) call cpsi_data_dealloc(next) 972 call cpsi_data_dealloc(psi2_tag) 973 call cpsi_data_dealloc(hml_tag) 974 call cpsi_data_dealloc(eig_tag) 975 call cpsi_data_dealloc(svec_tag) 976 value = BA_free_heap(dn(2)) 977 if (.not. value) 978 > call errquit('band_sd:freeing heap memory',0,MA_ERR) 979 980 981 call C3dB_pfft_end() 982 call cpsi_data_end() 983 call C3dB_end(1) 984 if ((control_gga().ge.10).and.(control_gga().le.200)) then 985 call mask_end() 986 call G_end() 987 call D3dB_end(1) 988 end if 989 990* |***************************| 991****************** report consumed cputime ********************** 992* |***************************| 993 if (taskid.eq.MASTER) then 994 CALL current_second(cpu4) 995 996 T1=CPU2-CPU1 997 T2=CPU3-CPU2 998 T3=CPU4-CPU3 999 T4=CPU4-CPU1 1000 AV=T2/dble(icount*it_in) 1001 write(6,*) 1002 write(6,*) '-----------------' 1003 write(6,*) 'cputime in seconds' 1004 write(6,*) 'prologue : ',T1 1005 write(6,*) 'main loop : ',T2 1006 write(6,*) 'epilogue : ',T3 1007 write(6,*) 'total : ',T4 1008 write(6,*) 'cputime/step: ',AV 1009 write(6,*) 1010 call nwpw_timing_print_final(.true.,(icount*it_in)) 1011 CALL nwpw_MESSAGE(4) 1012 end if 1013 1014 1015 call Parallel3d_Finalize() 1016 call Parallel_Finalize() 1017 band_sd = value 1018 return 1019 1020 1021*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 1022 1000 FORMAT(10X,'****************************************************') 1023 1010 FORMAT(10X,'* *') 1024 1020 FORMAT(10X,'* Car-Parrinello solid-state calculation *') 1025 1030 FORMAT(10X,'* [ steepest descent minimization ] *') 1026 1035 FORMAT(10x,'* [ NorthWest Chemistry implementation ] *') 1027 1040 FORMAT(10X,'* version #1.00 03/14/09 *') 1028 1041 FORMAT(10X,'* This code was developed by Eric J. Bylaska *') 1029 1042 FORMAT(10X,'* *') 1030 1043 FORMAT(10X,'* *') 1031 1100 FORMAT(//) 1032 1110 FORMAT(10X,'================ BAND input data ===================') 1033 1111 FORMAT(/' number of processors used:',I16) 1034 1112 FORMAT( ' parallel mapping : 1d slab') 1035 1113 FORMAT( ' parallel mapping : 2d hilbert') 1036 1114 FORMAT( ' parallel mapping : balanced') 1037 1115 FORMAT(/' options:') 1038 1116 FORMAT( ' parallel mapping : not balanced') 1039 1117 FORMAT( ' processor grid :',I4,' x',I4,' x',I4) 1040 1118 FORMAT( ' parallel mapping : 2d hcurve') 1041 1120 FORMAT(5X,' ionic motion = ',A) 1042 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 1043 1130 FORMAT(5X,' electron spin = ',A) 1044 1131 FORMAT(5X,' exchange-correlation = ',A) 1045 1132 FORMAT(5X,' using fractional occupation') 1046 1135 FORMAT(/' The masses of QM H atoms converted to 2.0 amu.', 1047 > /' To turn off this default', 1048 > ' set nwpw:makehmass2 .false.') 1049 1140 FORMAT(/' elements involved in the cluster:') 1050 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I1) 1051 1151 FORMAT(5X,' cutoff =',4F8.3) 1052 1152 FORMAT(12X,' highest angular component : ',i3) 1053 1153 FORMAT(12X,' local potential used : ',i3) 1054 1154 FORMAT(12X,' number of non-local projections: ',i3) 1055 1155 FORMAT(12X,' semicore corrections included : ', 1056 > F6.3,' (radius) ',F6.3,' (charge)') 1057 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 1058 1157 FORMAT(12X,' comment : ',A) 1059 1158 FORMAT(12X,' pseudpotential type : ',i3) 1060 1061 1159 FORMAT(/' total charge=',F8.3) 1062 1160 FORMAT(/' atomic composition:') 1063 1170 FORMAT(7(5X,A2,':',I3)) 1064 1180 FORMAT(/' initial position of ions:') 1065 1190 FORMAT(5X, I4, A5, ' (',3F11.5,' ) - atomic mass= ',F7.3,' ',A) 1066 1191 FORMAT(5X, I4, A5, ' (',3F11.5, 1067 > ' ) - atomic mass= ',F6.3,' - fixed ',A) 1068 1193 FORMAT(5X, I4, A5, ' (',3F11.5, 1069 > ' ) - atomic mass= ',F7.3,' - z fixed') 1070 1194 FORMAT(5X, I4, A5, ' (',3F11.5, 1071 > ' ) - atomic mass= ',F7.3,A) 1072 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 1073 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 1074 1211 FORMAT(5X,' number of constraints = ', I6,' ( DOF = ',I6,' )' ) 1075 1220 FORMAT(/' number of electrons: spin up=',F6.2, 16x, 1076 > ' down=',F6.2,A) 1077c 1220 FORMAT(/' number of electrons: spin up=',I6, 1078c > ' (',I4,' per task)', 1079c > ' down=',I6, 1080c > ' (',I4,' per task)', 1081c > A) 1082 1221 FORMAT( ' number of orbitals : spin up=',I6, 1083 > ' (',I4,' per task)', 1084 > ' down=',I6, 1085 > ' (',I4,' per task)', 1086 > A) 1087 1230 FORMAT(/' supercell:') 1088 1231 FORMAT(5x,' volume : ',F10.1) 1089 1232 FORMAT(/5x,' lattice: a=',f8.3,' b=',f8.3,' c=',f8.3, 1090 > /5x,' alpha=',f8.3,' beta=',f8.3,' gamma=',f8.3) 1091 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 1092 1242 FORMAT(5x,' a2=<',3f8.3,' >') 1093 1243 FORMAT(5x,' a3=<',3f8.3,' >') 1094 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >') 1095 1245 FORMAT(5x,' b2=<',3f8.3,' >') 1096 1246 FORMAT(5x,' b3=<',3f8.3,' >') 1097 1098 1249 FORMAT(/' computational grids:') 1099 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 1100 & '( ',I8,' waves ',I8,' per task)') 1101 1251 FORMAT(5X,' wavefnc ',I3,' cutoff=',F7.3, 1102 & ' fft=',I3,'x',I3,'x',I3, 1103 & '( ',I8,' waves ',I8,' per task)') 1104 1252 FORMAT(5x, ' wavefunction grids not printed - ', 1105 > 'number of k-point is very large') 1106 1107 1255 FORMAT(/' brillouin zone:') 1108 1256 FORMAT(5x,' number of zone points:',I3) 1109 1257 FORMAT(5x,' weight=',f8.3,' ks=<',3f8.3,' >, k=<',3f8.3,'>') 1110 1258 FORMAT(5x,' number of k-point is very large or distributed') 1111 1112 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,' and',I3) 1113 1261 FORMAT(5X,' madelung=',f14.8) 1114 1270 FORMAT(/' technical parameters:') 1115 1271 FORMAT(5x, ' translation constrained') 1116 1272 FORMAT(5x, ' rotation constrained') 1117 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1) 1118 1281 FORMAT(5X, ' maximum iterations =',I10, 1119 > ' ( ',I4,' inner ',I6,' outer )') 1120 1290 FORMAT(5X, ' tolerance=',E8.3,' (energy)',E12.3, 1121 & ' (electron)',E12.3,' (ion)') 1122 1297 FORMAT(/' fractional smearing parameters:') 1123 1298 FORMAT(5X, ' smearing algorithm = ',A) 1124 1299 FORMAT(5X, ' smearing parameter = ',E9.3,' (',F7.1,' K)') 1125 1300 FORMAT(//) 1126 1305 FORMAT(10X,'================ iteration =========================') 1127 1310 FORMAT(I8,E20.10,3E15.5) 1128 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 1129 1330 FORMAT(/' comparison between hamiltonian and lambda matrix') 1130 1331 FORMAT(/' Elements of Hamiltonian matrix (up/restricted)') 1131 1332 FORMAT(/' Elements of Hamiltonian matrix (down)') 1132 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7) 1133 1341 FORMAT(I3,2I3,' H=',E16.6) 1134 1350 FORMAT(/' orthonormality') 1135 1360 FORMAT(I3,2I3,'(',2E18.7,')') 1136 1370 FORMAT(I3) 1137 1380 FORMAT(' ''',a,'''',I4) 1138 1390 FORMAT(I3) 1139 1400 FORMAT(I3,3E18.8/3X,3E18.8) 1140 1410 FORMAT(10X,'============= summary of results =================') 1141 1420 FORMAT( ' final position of ions:') 1142 1430 FORMAT(//' total energy :',E19.10,' (',E15.5,'/ion)') 1143 1431 FORMAT(/' QM Energies') 1144 1432 FORMAT( '------------') 1145 1433 FORMAT( ' total QM energy :',E19.10,' (',E15.5,'/ion)') 1146 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)') 1147 1450 FORMAT( ' hartree energy :',E19.10,' (',E15.5,'/electron)') 1148 1455 FORMAT( ' SIC-hartree energy :',E19.10,' (',E15.5,'/electron)') 1149 1456 FORMAT( ' SIC-exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1150 1457 FORMAT( ' HF exchange energy :',E19.10,' (',E15.5,'/electron)') 1151 1460 FORMAT( ' exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1152 1470 FORMAT( ' ion-ion energy :',E19.10,' (',E15.5,'/ion)') 1153 1480 FORMAT(/' K.S. kinetic energy :',E19.10,' (',E15.5,'/electron)') 1154 1490 FORMAT( ' K.S. V_l energy :',E19.10,' (',E15.5,'/electron)') 1155 1495 FORMAT( ' K.S. V_nl energy :',E19.10,' (',E15.5,'/electron)') 1156 1496 FORMAT( ' K.S. V_Hart energy :',E19.10,' (',E15.5,'/electron)') 1157 1497 FORMAT( ' K.S. V_xc energy :',E19.10,' (',E15.5,'/electron)') 1158 1498 FORMAT( ' Virial Coefficient :',E19.10) 1159 1499 FORMAT( ' K.S. SIC-hartree energy :',E19.10, 1160 > ' (',E15.5,'/electron)') 1161 1501 FORMAT( ' K.S. SIC-exc-corr energy :',E19.10, 1162 > ' (',E15.5,'/electron)') 1163 1502 FORMAT( ' K.S. HFX energy :',E19.10, 1164 > ' (',E15.5,'/electron)') 1165 1500 FORMAT(/' orbital energies:') 1166 1507 FORMAT(/' Fermi energy =',2(E18.7,' (',F8.3,'eV)')) 1167 1508 FORMAT(/' Brillouin zone point: ',i3, 1168 > /' weight=',f10.6, 1169 > /' k =<',3f8.3,'> . <b1,b2,b3> ', 1170 > /' =<',3f8.3,'>') 1171 1510 FORMAT(4(E18.7,' (',F8.3,'eV) occ=',F5.3)) 1172 1511 FORMAT(33x,"Spin(Sz,Sy,Sz)") 1173 1512 FORMAT(E18.7,' (',F8.3,' eV) (|s| =',F6.3, 1174 > ', s = <',F7.3,',',F7.3,',',F7.3,'> ) occ=',F5.3) 1175 1176 1600 FORMAT(/' Total BAND energy :',E19.10) 1177 1178 1700 FORMAT(/' QM/MM-pol-vib/CAV Energies') 1179 1701 FORMAT( ' --------------------------') 1180 1702 FORMAT( ' LJ energy :',E19.10) 1181 1703 FORMAT( ' Residual Coulomb energy:',E19.10) 1182 1704 FORMAT( ' MM Vibration energy :',E19.10) 1183 1705 FORMAT( ' MM Vibration energy :',E19.10) 1184 1706 FORMAT( ' (QM+MM)/Cavity energy :',E19.10) 1185 1186 1720 FORMAT(/' Dispersion energy :',E19.10) 1187 1188 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 1189 1190 9000 if (taskid.eq.MASTER) write(6,9010) ierr 1191 call Parallel_Finalize() 1192 1193 band_sd = value 1194 return 1195 END 1196