1* 2* $Id$ 3* 4 5* **************************************************************** 6* * * 7* * wvfnc_new_lcao_density * 8* * * 9* **************************************************************** 10 logical function wvfnc_new_lcao_density(oprint_in, 11 > wavefunction_filename, 12 > version,ngrid,unita, 13 > ispin,ne,oddelcfill) 14 15 implicit none 16 logical oprint_in 17 character*50 wavefunction_filename 18 integer version 19 integer ngrid(3) 20 real*8 unita(3,3) 21 integer ispin,ne(2) 22 logical oddelcfill 23 24#include "bafdecls.fh" 25#include "errquit.fh" 26 27 28* **** local variables **** 29 integer MASTER,taskid 30 parameter (MASTER=0) 31 32 real*8 ALPHA,talpha 33 parameter (ALPHA=0.30d0) 34 35 logical value,value2,oprint,field_exist 36 integer ms,i,ia,icharge,ne_excited(2),nn,l,m,n,ii 37 integer dn(2),dnall(2),phi1(2),rho(2),vc(2),xcp(2),xce(2) 38 integer psi1(2),occ1(2),smearoccupation 39 integer ee(2),eigs(2),neq(2),mapping1d 40 integer dng(2),tmp1(2),vall(2),v_field(2),vl_lr(2),vl(2) 41 integer nx,ny,nz 42 integer nbasis,nbasis2,n2ft3d,npack0,npack1,it 43 real*8 rho_error,rho_error_old,total_dn,sum 44 real*8 Etotal,Eorbs,Ehart,Eexc,Eion,EV,de,olde,pxc 45 real*8 cpu1,cpu2,cpu3,cpu4,cpu5,t1,t2,t3,t4,av 46 real*8 scal1,scal2,dv 47 character*50 filename 48 49* **** external functions **** 50 logical psp_semicore,aorbs_init,aorbs_readall 51 character*4 ion_atom_qm,ion_atom 52 integer aorbs_norbs,aorbs_nbasis,ion_nkatm,psp_lmax,aorbs_katm 53 integer psp_locp,psp_lmmax,ion_natm,ion_nion,ion_nkatm_qm 54 integer control_version,aorbs_l,aorbs_m,aorbs_ii 55 integer ewald_ncut,ewald_nshl3d 56 real*8 psp_ncore,psp_rcore,psp_rlocal,psp_rc,psp_zv 57 real*8 ewald_rcut,ewald_e,ion_ion_e,coulomb_e 58 real*8 ewald_mandelung 59 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 60 real*8 lattice_unitg,control_tole,control_tolc 61 external psp_semicore,aorbs_init,aorbs_readall 62 external ion_atom_qm,ion_atom 63 external aorbs_norbs,aorbs_nbasis,ion_nkatm,psp_lmax,aorbs_katm 64 external psp_locp,psp_lmmax,ion_natm,ion_nion,ion_nkatm_qm 65 external control_version,aorbs_l,aorbs_m,aorbs_ii 66 external ewald_ncut,ewald_nshl3d 67 external psp_ncore,psp_rcore,psp_rlocal,psp_rc,psp_zv 68 external ewald_rcut,ewald_e,ion_ion_e,coulomb_e 69 external ewald_mandelung 70 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 71 external lattice_unitg,control_tole,control_tolc 72 integer pack_nwave_all 73 integer control_gga 74 integer control_ngrid,pack_nwave 75 external pack_nwave_all 76 external control_gga 77 external control_ngrid,pack_nwave 78 character*12 control_boundry 79 external control_boundry 80 81 integer control_excited_ne,control_mapping1d 82 character*50 control_input_epsi 83 external control_excited_ne,control_mapping1d 84 external control_input_epsi 85 logical ion_chargeexist,ion_mmexist 86 external ion_chargeexist,ion_mmexist 87 88 89 90 91 Etotal = 0.0d0 92 oprint = oprint_in 93 94 call Parallel_taskid(taskid) 95 96 if ((taskid.eq.MASTER).and.oprint) call current_second(cpu1) 97 98 call D3dB_n2ft3d(1,n2ft3d) 99 call Pack_npack(0,npack0) 100 call Pack_npack(1,npack1) 101 102 field_exist = ion_chargeexist().or.ion_mmexist() 103 104 mapping1d = control_mapping1d() 105 call Dne_init(ispin,ne,mapping1d) 106 call Dneall_neq(neq) 107 108 call D3dB_nx(1,nx) 109 call D3dB_ny(1,ny) 110 call D3dB_nz(1,nz) 111 scal1 = 1.0d0/dble(nx*ny*nz) 112 113* **** initialize atomic orbitals **** 114 value = aorbs_init() 115 value = value.and.aorbs_readall() 116 if (.not.value) go to 101 117 nbasis = aorbs_nbasis() 118 nbasis2 = nbasis**2 119 120* **** basis set not big enough to meeting electronic filling **** 121 if ((nbasis.lt.ne(1)).or.(nbasis.lt.ne(2))) then 122 value = .false. 123 go to 101 124 end if 125 126 127* ***** allocate memory from heap memory **** 128 value = value.and. 129 > BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 130 > 'psi1',psi1(2),psi1(1)) 131 value = value.and. 132 > BA_alloc_get(mt_dbl,(2*n2ft3d),'dn',dn(2),dn(1)) 133 value = value.and. 134 > BA_alloc_get(mt_dbl,(2*n2ft3d),'dnall',dnall(2),dnall(1)) 135 value = value.and. 136 > BA_alloc_get(mt_dcpl,(npack0),'dng',dng(2),dng(1)) 137 value = value.and. 138 > BA_alloc_get(mt_dcpl,(npack1),'phi1',phi1(2),phi1(1)) 139 value = value.and. 140 > BA_alloc_get(mt_dbl,(2*nbasis),'ee',ee(2),ee(1)) 141 value = value.and. 142 > BA_alloc_get(mt_dbl,(n2ft3d),'rho',rho(2),rho(1)) 143 value = value.and. 144 > BA_alloc_get(mt_dbl,(n2ft3d),'vc',vc(2),vc(1)) 145 value = value.and. 146 > BA_alloc_get(mt_dbl,(2*n2ft3d),'xce',xce(2),xce(1)) 147 value = value.and. 148 > BA_alloc_get(mt_dbl,(2*n2ft3d),'xcp',xcp(2),xcp(1)) 149 value = value.and. 150 > BA_alloc_get(mt_dbl,(2*n2ft3d),'vall',vall(2),vall(1)) 151 value = value.and. 152 > BA_alloc_get(mt_dcpl,npack0, 153 > 'vl2',vl(2),vl(1)) 154 value = value.and. 155 > BA_alloc_get(mt_dbl,n2ft3d, 156 > 'vl_lr',vl_lr(2),vl_lr(1)) 157 value = value.and. 158 > BA_alloc_get(mt_dbl,n2ft3d, 159 > 'v_field',v_field(2),v_field(1)) 160 eigs(1) = ee(1) 161 eigs(2) = ee(1) + nbasis 162 163* ***** allocate matrices using GA **** 164 if (.not. value) then 165 go to 100 166 end if 167 168 if ((taskid.eq.MASTER).and.(oprint)) then 169 write(*,1000) 170 write(*,1010) 171 write(*,1020) 172 write(*,1010) 173 write(*,1030) 174 write(*,1010) 175 write(*,1040) 176 write(*,1010) 177 write(*,1000) 178 write(*,1110) 179 call nwpw_message(1) 180 181 182 write(6,1121) control_boundry(),control_version() 183 if (ispin.eq.1) write(6,1130) 'restricted' 184 if (ispin.eq.2) write(6,1130) 'unrestricted' 185 186 IF (control_gga().eq.0) THEN 187 write(6,1131) 'Vosko et al parameterization' 188 ELSE IF (control_gga().eq.10) THEN 189 write(6,1131) 190 > 'PBE96 (White and Bird) parameterization' 191 ELSE IF (control_gga().eq.11) THEN 192 write(6,1131) 193 > 'BLYP (White and Bird) parameterization' 194 ELSE IF (control_gga().eq.12) THEN 195 write(6,1131) 196 > 'revPBE (White and Bird) parameterization' 197 ELSE 198 write(6,1131) 'unknown parameterization' 199 call errquit('bad exchange_correlation',0, INPUT_ERR) 200 END IF 201 202 write(*,1117) 203 do ia=1,ion_nkatm_qm() 204 write(*,1118) ion_atom_qm(ia),aorbs_norbs(ia) 205 end do 206 write(*,1119) nbasis 207 208 do n = 1,nbasis 209 ii = aorbs_ii(n) 210 ia = aorbs_katm(n) 211 l = aorbs_l(ia,n) 212 m = aorbs_m(ia,n) 213 write(*,*) "basis function=",n," ii,ia,l,m=",ii,ia,l,m 214 end do 215 216 217 write(6,1140) 218 do ia = 1,ion_nkatm() 219 write(6,1150) ia,ion_atom(ia), 220 > psp_zv(ia),psp_lmax(ia) 221 write(6,1152) psp_lmax(ia) 222 write(6,1153) psp_locp(ia) 223 write(6,1154) psp_lmmax(ia) 224 if (control_version().eq.4) write(6,1156) psp_rlocal(ia) 225 if (psp_semicore(ia)) 226 > write(6,1155) psp_rcore(ia),psp_ncore(ia) 227 write(6,1151) (psp_rc(i,ia),i=0,psp_lmax(ia)) 228 end do 229 230 icharge = -(ne(1)+ne(ispin)) 231 do ia=1,ion_nkatm() 232 icharge = icharge + ion_natm(ia)*psp_zv(ia) 233 end do 234 write(6,1159) icharge 235 236 write(*,1220) ispin, ne(1), ne(ispin) 237 write(6,1230) 238 write(6,1241) lattice_unita(1,1), 239 > lattice_unita(2,1), 240 > lattice_unita(3,1) 241 write(6,1242) lattice_unita(1,2), 242 > lattice_unita(2,2), 243 > lattice_unita(3,2) 244 write(6,1243) lattice_unita(1,3), 245 > lattice_unita(2,3), 246 > lattice_unita(3,3) 247 write(6,1244) lattice_unitg(1,1), 248 > lattice_unitg(2,1), 249 > lattice_unitg(3,1) 250 write(6,1245) lattice_unitg(1,2), 251 > lattice_unitg(2,2), 252 > lattice_unitg(3,2) 253 write(6,1246) lattice_unitg(1,3), 254 > lattice_unitg(2,3), 255 > lattice_unitg(3,3) 256 write(6,1231) lattice_omega() 257 write(6,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 258 > pack_nwave_all(0),pack_nwave(0) 259 write(6,1251) lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 260 > pack_nwave_all(1),pack_nwave(1) 261 if (control_version().eq.3) then 262 write(6,1260) ewald_rcut(),ewald_ncut() 263 write(6,1261) ewald_mandelung() 264 end if 265 write(6,1270) 266 write(6,1300) 267 write(6,1305) 268 call util_flush(6) 269 270 write(6,*) 271 end if 272 if ((taskid.eq.MASTER).and.oprint) call current_second(cpu2) 273 274 275 276* ******************************** 277* **** generate phaze factors **** 278* ******************************** 279 call phafac() 280 if (psp_semicore(0)) call semicore_density_update() 281 282 283* ******************************** 284* **** intialize the density **** 285* ******************************** 286 call lcao_init_dn(ispin,ne,n2ft3d,dbl_mb(dn(1)),dcpl_mb(phi1(1))) 287 288 289* ********************** 290* **** generate dng **** 291* ********************** 292 value = BA_push_get(mt_dbl,n2ft3d,'tmp1',tmp1(2),tmp1(1)) 293 if (.not. value) call errquit( 294 > 'electron_gen_dng_dnall: out of stack memory',0, MA_ERR) 295 296 297 call D3dB_rr_Sum(1,dbl_mb(dn(1)), 298 > dbl_mb(dn(1)+(ispin-1)*n2ft3d), 299 > dbl_mb(tmp1(1))) 300 call D3dB_r_SMul1(1,scal1,dbl_mb(tmp1(1))) 301 call D3dB_rc_fft3f(1,dbl_mb(tmp1(1))) 302 call Pack_c_pack(0,dbl_mb(tmp1(1))) 303 call Pack_c_Copy(0,dbl_mb(tmp1(1)),dcpl_mb(dng(1))) 304 305* ******************************************************** 306* **** generate dnall - used for semicore corrections **** 307* ******************************************************** 308 if (psp_semicore(0)) then 309 call semicore_density(dbl_mb(tmp1(1))) 310 call D3dB_r_SMul1(1,0.5d0,dbl_mb(tmp1(1))) 311 else 312 call D3dB_r_Zero(1,dbl_mb(tmp1(1))) 313 end if 314 do ms=1,ispin 315 call D3dB_rr_Sum(1,dbl_mb(dn(1)+(ms-1)*n2ft3d), 316 > dbl_mb(tmp1(1)), 317 > dbl_mb(dnall(1)+(ms-1)*n2ft3d)) 318 end do 319 320* ***************************** 321* **** generate potentials **** 322* ***************************** 323 if (control_version().eq.3) then 324 call coulomb_v(dcpl_mb(dng(1)),dbl_mb(vc(1))) 325 end if 326 327 if (control_version().eq.4) then 328 call D3dB_rr_Sum(1,dbl_mb(dn(1)), 329 > dbl_mb(dn(1)+(ispin-1)*n2ft3d), 330 > dbl_mb(tmp1(1))) 331 332 call coulomb2_v(dbl_mb(tmp1(1)),dbl_mb(vc(1))) 333 end if 334 335* **** xc potential **** 336 call v_bwexc_all(control_gga(),n2ft3d,ispin,dbl_mb(dnall(1)), 337 > dbl_mb(xcp(1)),dbl_mb(xce(1))) 338 value = BA_pop_stack(tmp1(2)) 339 if (.not.value) call errquit( 340 > 'electron_gen_dng_dnall: poping stack',1, MA_ERR) 341 342 343 344 scal2 = 1.0d0/lattice_omega() 345 346 if (control_version().eq.3) then 347 348* **** add up k-space potentials, vall = scal2*vl + vc **** 349 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 350 > dbl_mb(vall(1))) 351 352 call Pack_cc_Sum2(0,dcpl_mb(vc(1)),dbl_mb(vall(1))) 353 354* **** fourier transform k-space potentials **** 355 call Pack_c_unpack(0,dbl_mb(vall(1))) 356 call D3dB_cr_fft3b(1,dbl_mb(vall(1))) 357 if (field_exist) 358 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 359 360 else 361 362* **** add up k-space potentials, vall = scal2*vsr_l **** 363 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 364 > dbl_mb(vall(1))) 365 366* **** fourier transform k-space potentials **** 367 call Pack_c_unpack(0,dbl_mb(vall(1))) 368 call D3dB_cr_fft3b(1,dbl_mb(vall(1))) 369 370 call D3dB_rr_Sum2(1,dbl_mb(vl_lr(1)),dbl_mb(vall(1))) 371 call D3dB_rr_Sum2(1,dcpl_mb(vc(1)),dbl_mb(vall(1))) 372 if (field_exist) 373 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 374 375 end if 376 377 if (ispin.eq.2) then 378 call D3dB_rr_Sum(1,dbl_mb(vall(1)), 379 > dbl_mb(xcp(1) +n2ft3d), 380 > dbl_mb(vall(1)+n2ft3d)) 381 call D3dB_r_Zero_Ends(1,dbl_mb(vall(1)+n2ft3d)) 382 end if 383 call D3dB_rr_Sum2(1,dbl_mb(xcp(1)),dbl_mb(vall(1))) 384 call D3dB_r_Zero_Ends(1,dbl_mb(vall(1))) 385 386 387* ******************************* 388* *** read psi1 wavefunctions *** 389* ******************************* 390 write(*,*) "HERA" 391 occ1(1) = rho(1) 392 call psi_get_ne_occupation(ispin,ne,smearoccupation) 393 if (smearoccupation.gt.0) then 394 write(*,*) "HERAaa" 395 value = value.and. 396 > BA_alloc_get(mt_dbl,(ne(1)+ne(2)),'occ1',occ1(2),occ1(1)) 397 end if 398 write(*,*) "HERB,ispin,ne,smearocc=",ispin,ne,smearoccupation 399 call psi_read(ispin,ne,dcpl_mb(psi1(1)), 400 > smearoccupation,dbl_mb(occ1(1))) 401 402 write(*,*) "HERC,ispin,ne,smearocc=",ispin,ne,smearoccupation 403* ******************************** 404* *** write psi1 wavefunctions *** 405* ******************************** 406 call psi_write(ispin,ne,dcpl_mb(psi1(1)), 407 > smearoccupation,dbl_mb(occ1(1))) 408 409 410 411 412 413 414 rho_error = 10000.0d0 415 if ((taskid.eq.MASTER).and.oprint) call current_second(cpu3) 416 if ((taskid.eq.MASTER).and.oprint) CALL nwpw_MESSAGE(8) 417 418 419 50 if ((taskid.eq.MASTER).and.oprint) call current_second(cpu4) 420 if ((taskid.eq.MASTER).and.oprint) CALL nwpw_MESSAGE(3) 421 if ((taskid.eq.MASTER).and.oprint) write(*,*) 422 if ((taskid.eq.MASTER).and.oprint) write(*,*) 423 424 425c 50 call lcao_diis_end() 426c call lcao_write_psi(wavefunction_filename, 427c > version, 428c > ngrid, 429c > unita, 430c > ispin,ne, 431c > psimatrix,dcpl_mb(phi2(1))) 432c 433 434* **** write out excited orbitals **** 435 ne_excited(1) = 0 436 ne_excited(2) = 0 437 ne_excited(1) = control_excited_ne(1) 438 if (ispin.eq.2) ne_excited(2) = control_excited_ne(2) 439 filename = control_input_epsi() 440 441c if ((ne_excited(1)+ne_excited(2)).gt.0) 442c > call lcao_write_epsi(filename, 443c > version, 444c > ngrid, 445c > unita, 446c > ispin,ne,ne_excited, 447c > psimatrix,dcpl_mb(phi2(1))) 448c 449 450*::::::::::::::::: report summary of results ::::::::::::::::::::::: 451 452 if ((taskid.eq.MASTER).and.oprint) then 453 write(*,*) 454 write(*,1410) 455 end if 456 457C **** caluclate number of Density **** 458 do ms=1,ispin 459 call D3dB_r_dsum(1,dbl_mb(dn(1)+(ms-1)*n2ft3d),total_dn) 460 total_dn = total_dn*lattice_omega() 461 total_dn = total_dn/dble(ngrid(1)*ngrid(2)*ngrid(3)) 462 if ((taskid.eq.MASTER).and.oprint) then 463 write(*,*) 'Spin=',MS,' Total Orbital Density:', total_dn 464 end if 465 end do 466 467 468 469 470 471*::::::::::::::::::: report consumed cputime ::::::::::::::::::::::: 472 if ((taskid.eq.MASTER).and.oprint) call current_second(cpu5) 473 474 475 if ((taskid.eq.MASTER).and.oprint) then 476 t1 = cpu2 - cpu1 477 t2 = cpu3 - cpu2 478 t3 = cpu4 - cpu3 479 t4 = cpu5 - cpu1 480 WRITE(*,*) 481 WRITE(*,*) '----------------------' 482 WRITE(*,*) 'Proglogue : ',t1 483 WRITE(*,*) 'Setup Matrices : ',t2 484 WRITE(*,*) 'Convergence Loop : ',t3 485 WRITE(*,*) 'Total : ',t4 486 WRITE(*,*) 'CpuTime/Step : ',av 487 end if 488 if ((taskid.eq.MASTER).and.oprint) CALL nwpw_MESSAGE(4) 489 490 491 492 493 494* ***** free heap memory **** 495 100 continue 496 value2 = BA_free_heap(vc(2)) 497 value2 = value2.and.BA_free_heap(xce(2)) 498 value2 = value2.and.BA_free_heap(xcp(2)) 499 value2 = value2.and.BA_free_heap(vall(2)) 500 value2 = value2.and.BA_free_heap(vl_lr(2)) 501 value2 = value2.and.BA_free_heap(vl(2)) 502 value2 = value2.and.BA_free_heap(v_field(2)) 503 value2 = value2.and.BA_free_heap(rho(2)) 504 value2 = value2.and.BA_free_heap(dn(2)) 505 value2 = value2.and.BA_free_heap(psi1(2)) 506 value2 = value2.and.BA_free_heap(dng(2)) 507 value2 = value2.and.BA_free_heap(phi1(2)) 508 value2 = value2.and.BA_free_heap(ee(2)) 509 call Dne_end() 510 101 call aorbs_end() 511 512 value = .false. 513 wvfnc_new_lcao_density = value 514 return 515 516 1000 FORMAT(10X,'****************************************************') 517 1010 FORMAT(10X,'* *') 518 1020 FORMAT(10X,'* LCAO Calculations *') 519 1030 FORMAT(10X,'* [ NorthWest Chemistry implementation ] *') 520 1040 FORMAT(10X,'* Version #2.00 *') 521 1100 FORMAT(//) 522 1110 FORMAT(10X,'================ input data ========================') 523 1111 FORMAT(' Reading in ELCIN file(1-yes, 0-no): ',I1) 524 1115 FORMAT(/' Pseudopotentials used:') 525 1116 FORMAT(5X,A,' # of L: ',I1) 526 1117 FORMAT(/' Atomic orbitals used:') 527 1118 FORMAT(5X,A,' # of orbitals: ',I6) 528 1119 FORMAT(/' Number of basis functions:',I4) 529 1121 FORMAT(/5X,' boundry conditions = ',A,'(version', I1,')') 530 1130 FORMAT(5X,' electron spin = ',A) 531 1131 FORMAT(5X,' exchange-correlation = ',A) 532 1140 FORMAT(/' elements involved in the cluster:') 533 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I1) 534 1151 FORMAT(5X,' cutoff =',4F8.3) 535 1152 FORMAT(12X,' highest angular component : ',i2) 536 1153 FORMAT(12X,' local potential used : ',i2) 537 1154 FORMAT(12X,' number of non-local projections: ',i2) 538 1155 FORMAT(12X,' semicore corrections included : ', 539 > F6.3,' (radius) ',F6.3,' (charge)') 540 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 541 1159 FORMAT(/' total charge=',I2) 542 543 1160 FORMAT(/' ATOMIC COMPOSITION:') 544 1170 FORMAT(7(5X,A2,':',I3)) 545 1180 FORMAT(/' INITIAL POSITION OF IONS:') 546 1190 FORMAT(5X, I4, A4 ,' (',3F11.5,' )') 547 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 548 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 549 1220 FORMAT(' # of Electrons(ispin=',I1,') :',I4,' up',I4,' down') 550 1230 FORMAT(/' supercell:') 551 1231 FORMAT(5x,' volume : ',F10.1) 552 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 553 1242 FORMAT(5x,' a2=<',3f8.3,' >') 554 1243 FORMAT(5x,' a3=<',3f8.3,' >') 555 1244 FORMAT(5x,' b1=<',3f8.3,' >') 556 1245 FORMAT(5x,' b2=<',3f8.3,' >') 557 1246 FORMAT(5x,' b3=<',3f8.3,' >') 558 559 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 560 & '( ',I8,' waves ',I8,' per task)') 561 1251 FORMAT(5X,' wavefnc cutoff=',F7.3,' fft=',I3,'x',I3,'x',I3, 562 & '( ',I8,' waves ',I8,' per task)') 563 564 1260 FORMAT(5X,' ewald summation: cut radius=',F8.2,' and',I3) 565 1261 FORMAT(5X,' madelung=',f14.8) 566 1270 FORMAT(/' technical parameters:') 567 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictacious mass=',F10.1) 568 1290 FORMAT(5X, ' tolerance=',E8.3,' (energy)',E12.3, 569 & ' (electron)',E12.3,' (ion)') 570 1300 FORMAT(//) 571 572 1305 FORMAT(10X,'================ ITERATION =========================') 573 1310 FORMAT(I8,E20.10,2E15.5) 574 1320 FORMAT(' NUMBER OF ELECTRONS: SPIN UP=',F11.5,' DOWN=',F11.5) 575 1330 FORMAT(/' COMPARISON BETWEEN HAMILTONIAN AND LAMBDA MATRIX') 576 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7) 577 1350 FORMAT(/' ORTHONORNALITY') 578 1360 FORMAT(I3,2I3,E18.7) 579 1370 FORMAT(I3) 580 1380 FORMAT(' ''',a,'''',I4) 581 1390 FORMAT(I3) 582 1400 FORMAT(I3,3E18.8/3X,3E18.8) 583 1410 FORMAT(10X,'============= SUMMARY OF RESULTS =================') 584 1420 FORMAT( ' FINAL POSITION OF IONS:') 585 1430 FORMAT(/' TOTAL ENERGY :',E19.10,' (',E15.5,'/ION)') 586 1440 FORMAT( ' TOTAL ORBITAL ENERGY:',E19.10,' (',E15.5,'/ION)') 587 1450 FORMAT( ' HARTREE ENERGY :',E19.10,' (',E15.5,'/ELECTRON)') 588 1460 FORMAT( ' EXC-CORR ENERGY :',E19.10,' (',E15.5,'/ELECTRON)') 589 1470 FORMAT( ' ION-ION ENERGY :',E19.10,' (',E15.5,'/ION)') 590 1480 FORMAT( ' PXC ENERGY :',E19.10,' (',E15.5,'/ION)') 591 1500 FORMAT(/' ORBITAL ENERGIES:') 592 1510 FORMAT(2(E18.7,' (',F8.3,'eV)')) 593 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 594 595 end 596 597 598 599