1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! This module contains the following active subroutines: 3! SCFatom_Init,SCFatom,Orbit_Init,NC_Init,SC_Init,FC_Init,Potential_Init, 4! Set_Valence 5! The following subroutines need revision to be useful: 6! dump_aeatom,load_aeatom 7!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 8 9#if defined HAVE_CONFIG_H 10#include "config.h" 11#endif 12 13MODULE AEatom 14 15 USE io_tools 16 USE atomdata 17 USE excor 18 USE exx_mod 19 USE general_mod 20 USE globalmath 21 USE gridmod 22 USE hf_mod 23 USE ldagga_mod 24 USE input_dataset_mod 25 26 IMPLICIT NONE 27 28 REAL(8), PARAMETER, PRIVATE :: linh=0.0025d0,logh=0.020d0,lor00=1.d-5 29 REAL(8), PARAMETER, PRIVATE :: small0=1.d-13,rimix=0.5d0,worst=1.d-8 30 INTEGER, PARAMETER, PRIVATE :: niter=1000,mxloop=500 31 REAL(8), PARAMETER, PRIVATE :: conv1=4.d13,conv2=3.d13,conv3=2.d13,conv4=1.d13 32 REAL(8), PRIVATE :: electrons 33 34 TYPE (PotentialInfo) ,TARGET :: AEPot 35 TYPE (OrbitInfo) ,TARGET :: AEOrbit 36 TYPE (SCFInfo) ,TARGET :: AESCF 37 38 TYPE (PotentialInfo) ,TARGET :: FCPot 39 TYPE (OrbitInfo) ,TARGET :: FCOrbit 40 TYPE (SCFInfo) ,TARGET :: FCSCF 41 TYPE(FCInfo),TARGET :: FC 42 43 TYPE (Gridinfo) ,TARGET :: Grid 44 45 TYPE (PotentialInfo) ,PRIVATE, POINTER :: PotPtr 46 TYPE (OrbitInfo) ,PRIVATE, POINTER :: OrbitPtr 47 TYPE (SCFInfo) ,PRIVATE, POINTER :: SCFPtr 48 49CONTAINS 50 51!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 52!! SCFatom_Init !!!! 53!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 54 55 SUBROUTINE SCFatom_Init() 56 57 ! program to calculate the self-consistent density functional 58 ! atom ground state for atom with atomic number nz 59 ! for self-consistent potential rv 60 61 REAL(8) :: xocc,qf,small,zeff,hval 62 REAL(8) :: qcal,rescale,nzeff,h,r0 63 INTEGER :: icount,i,j,k,it,start,np,ierr 64 INTEGER :: is,ip,id,jf,ig,io,l,nfix,ir,kappa 65 INTEGER :: nps,npp,npd,npf,npg 66 INTEGER, ALLOCATABLE :: nl(:,:) 67 68! Various initializations 69 AEPot%v0=0.d0;AEPot%v0p=0.d0;AEPot%Nv0=0;AEPot%Nv0p=0 70 frozencorecalculation=.FALSE.;setupfrozencore=.false. 71 gaussianshapefunction=.FALSE.;besselshapefunction=.FALSE. 72 73! Atomic symbol and atomic number (from input dataset) 74 AEPot%sym=input_dataset%atomic_symbol 75 AEPot%nz=input_dataset%atomic_charge 76 AEPot%zz=AEPot%nz 77 78! Relativistic, point-nucleus, HF data (from input dataset) 79 scalarrelativistic=input_dataset%scalarrelativistic 80 diracrelativistic=input_dataset%diracrelativistic 81 HFpostprocess=input_dataset%HFpostprocess 82 finitenucleus=input_dataset%finitenucleus 83 AEPot%finitenucleusmodel=input_dataset%finitenucleusmodel 84 85! Grid data (from input dataset) 86 IF (TRIM(input_dataset%gridkey)=='LINEAR') THEN 87 hval=input_dataset%gridmatch/(input_dataset%gridpoints-1) 88 CALL InitGrid(Grid,hval,input_dataset%gridrange) 89 ELSEIF (TRIM(input_dataset%gridkey)=='LOGGRID') THEN 90 hval=logh 91 CALL findh(AEPot%zz,input_dataset%gridmatch,input_dataset%gridpoints,hval,r0) 92 CALL InitGrid(Grid,hval,input_dataset%gridrange,r0=r0) 93 ELSEIF (TRIM(input_dataset%gridkey)=='LOGGRID4') THEN 94 hval=logh 95 CALL findh_given_r0(AEPot%zz,input_dataset%gridmatch,lor00,& 96& input_dataset%gridpoints,hval) 97 CALL InitGrid(Grid,hval,input_dataset%gridrange,r0=lor00/AEPot%zz) 98 ENDIF 99 100 CALL InitPot(AEPot,Grid%n) 101 CALL Get_Nuclearpotential(Grid,AEPot) 102 103! Logderiv data 104 minlogderiv=input_dataset%minlogderiv 105 maxlogderiv=input_dataset%maxlogderiv 106 nlogderiv=input_dataset%nlogderiv 107 108! Bound state solver 109 BDsolve=input_dataset%BDsolve 110 111! Exchange-correlation/HF keyword (from input dataset) 112 exctype=input_dataset%exctype 113 localizedcoreexchange=input_dataset%localizedcoreexchange 114 ColleSalvetti=.FALSE. 115 IF (TRIM(input_dataset%exctype)=='EXX'.or.& 116& TRIM(input_dataset%exctype)=='EXXKLI'.or.& 117& TRIM(input_dataset%exctype)=='EXXOCC') THEN 118 CALL EXX_Input_Settings(input_dataset%fixed_zero,input_dataset%fixed_zero_index) 119 ELSE IF (TRIM(input_dataset%exctype)=='EXXCS') THEN 120 CALL EXX_Input_Settings(input_dataset%fixed_zero,input_dataset%fixed_zero_index) 121 ColleSalvetti=.TRUE. 122 ELSE IF (TRIM(input_dataset%exctype)=='HF'.or.& 123& TRIM(input_dataset%exctype)=='HFV') THEN 124 ELSE 125 CALL initexch 126 ENDIF 127 128! Electronic configuration of atom 129 WRITE(STD_OUT,'(a,f6.2)') ' Calculation for atomic number = ',AEPot%zz 130 call flush_unit(std_out) 131 132 nps=input_dataset%np(1) 133 npp=input_dataset%np(2) 134 npd=input_dataset%np(3) 135 npf=input_dataset%np(4) 136 npg=input_dataset%np(5) 137 WRITE(STD_OUT,'(5i4)') nps,npp,npd,npf,npg 138 139 i=MAX(nps,npp,npd,npf,npg) 140 j=nps 141 IF(npp>0) j=j+npp-1 142 IF(npd>0) j=j+npd-2 143 IF(npf>0) j=j+npf-3 144 IF(npg>0) j=j+npg-4 145 146 IF (diracrelativistic) j=j+j ! need more orbitals 147 CALL InitOrbit(AEOrbit,j,Grid%n,exctype) 148 AEOrbit%nps=nps;AEOrbit%npp=npp;AEOrbit%npd=npd 149 AEOrbit%npf=npf;AEOrbit%npg=npg 150 ALLOCATE(nl(i,j));nl=0 151 152 icount=0 153 IF (.not.diracrelativistic) THEN 154 IF (nps.GT.0) THEN 155 DO is=1,nps 156 icount=icount+1 157 nl(is,1)=icount 158 AEOrbit%occ(icount)=2.d0 159 AEOrbit%np(icount)=is 160 AEOrbit%l(icount)=0 161 ENDDO 162 ENDIF 163 IF (npp.GT.1) THEN 164 DO ip=2,npp 165 icount=icount+1 166 nl(ip,2)=icount 167 AEOrbit%occ(icount)=6.d0 168 AEOrbit%np(icount)=ip 169 AEOrbit%l(icount)=1 170 ENDDO 171 ENDIF 172 IF (npd.GT.2) THEN 173 DO id=3,npd 174 icount=icount+1 175 nl(id,3)=icount 176 AEOrbit%occ(icount)=10.d0 177 AEOrbit%np(icount)=id 178 AEOrbit%l(icount)=2 179 ENDDO 180 ENDIF 181 IF (npf.GT.3) THEN 182 DO jf=4,npf 183 icount=icount+1 184 nl(jf,4)=icount 185 AEOrbit%occ(icount)=14.d0 186 AEOrbit%np(icount)=jf 187 AEOrbit%l(icount)=3 188 ENDDO 189 ENDIF 190 IF(npg.GT.4) THEN 191 DO ig=5,npg 192 icount=icount+1 193 nl(ig,5)=icount 194 AEOrbit%occ(icount)=18.d0 195 AEOrbit%np(icount)=ig 196 AEOrbit%l(icount)=4 197 ENDDO 198 ENDIF 199 AEOrbit%norbit=icount 200 AEOrbit%nps=nps 201 AEOrbit%npp=npp 202 AEOrbit%npd=npd 203 AEOrbit%npf=npf 204 AEOrbit%npg=npg 205 IF (AEOrbit%norbit/=input_dataset%norbit) STOP 'Inconsistent number of orbitals!' 206 WRITE(STD_OUT,*) AEOrbit%norbit, ' orbitals will be calculated' 207 208 WRITE(STD_OUT,*)' Below are listed the default occupations ' 209 WRITE(STD_OUT,"(' n l occupancy')") 210 DO io=1,AEOrbit%norbit 211 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,1e15.7)') & 212 & AEOrbit%np(io),AEOrbit%l(io),AEOrbit%occ(io) 213 ENDDO 214 215 !Corrected occupations (from input dataset) 216 DO io=1,input_dataset%norbit_mod 217 l=input_dataset%orbit_mod_l(io) 218 ip=input_dataset%orbit_mod_n(io) 219 xocc=input_dataset%orbit_mod_occ(io) 220 nfix=nl(ip,l+1) 221 IF (nfix<=0.OR.nfix>AEOrbit%norbit) THEN 222 WRITE(STD_OUT,*) 'error in occupations -- ip,l,xocc: ',ip,l,xocc,nfix,AEOrbit%norbit 223 STOP 224 ENDIF 225 AEOrbit%occ(nfix)=xocc 226 END DO 227 228 WRITE(STD_OUT,*) ' Corrected occupations are: ' 229 WRITE(STD_OUT,"(' n l occupancy')") 230 electrons=0.d0 231 DO io=1,AEOrbit%norbit 232 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,1e15.7)') & 233 & AEOrbit%np(io),AEOrbit%l(io),AEOrbit%occ(io) 234 electrons=electrons+AEOrbit%occ(io) 235 ENDDO 236 ENDIF ! scalarrelativistic 237 238 IF (diracrelativistic) THEN 239 icount=0 240 deallocate(nl) 241 i=MAX(nps,npp,npd,npf,npg) 242 allocate(nl(i,-5:5)) 243 nl=0 244 IF (nps.GT.0) THEN 245 DO is=1,nps 246 icount=icount+1 247 nl(is,-1)=icount 248 AEOrbit%occ(icount)=2.d0 249 AEOrbit%np(icount)=is 250 AEOrbit%l(icount)=0 251 AEOrbit%kappa(icount)=-1 252 ENDDO 253 ENDIF 254 IF (npp.GT.1) THEN 255 DO ip=2,npp 256 icount=icount+1 257 nl(ip,1)=icount 258 AEOrbit%occ(icount)=2.d0 259 AEOrbit%np(icount)=ip 260 AEOrbit%l(icount)=1 261 AEOrbit%kappa(icount)=1 262 ENDDO 263 DO ip=2,npp 264 icount=icount+1 265 nl(ip,-2)=icount 266 AEOrbit%occ(icount)=4.d0 267 AEOrbit%np(icount)=ip 268 AEOrbit%l(icount)=1 269 AEOrbit%kappa(icount)=-2 270 ENDDO 271 ENDIF 272 IF (npd.GT.2) THEN 273 DO id=3,npd 274 icount=icount+1 275 nl(id,2)=icount 276 AEOrbit%occ(icount)=4.d0 277 AEOrbit%np(icount)=id 278 AEOrbit%l(icount)=2 279 AEOrbit%kappa(icount)=2 280 ENDDO 281 DO id=3,npd 282 icount=icount+1 283 nl(id,-3)=icount 284 AEOrbit%occ(icount)=6.d0 285 AEOrbit%np(icount)=id 286 AEOrbit%l(icount)=2 287 AEOrbit%kappa(icount)=-3 288 ENDDO 289 ENDIF 290 IF (npf.GT.3) THEN 291 DO jf=4,npf 292 icount=icount+1 293 nl(jf,3)=icount 294 AEOrbit%occ(icount)=6.d0 295 AEOrbit%np(icount)=jf 296 AEOrbit%l(icount)=3 297 AEOrbit%kappa(icount)=3 298 ENDDO 299 DO jf=4,npf 300 icount=icount+1 301 nl(jf,-4)=icount 302 AEOrbit%occ(icount)=8.d0 303 AEOrbit%np(icount)=jf 304 AEOrbit%l(icount)=3 305 AEOrbit%kappa(icount)=-4 306 ENDDO 307 ENDIF 308 IF(npg.GT.4) THEN 309 DO ig=5,npg 310 icount=icount+1 311 nl(ig,4)=icount 312 AEOrbit%occ(icount)=8.d0 313 AEOrbit%np(icount)=ig 314 AEOrbit%l(icount)=4 315 AEOrbit%kappa(icount)=4 316 ENDDO 317 DO ig=5,npg 318 icount=icount+1 319 nl(ig,-5)=icount 320 AEOrbit%occ(icount)=10.d0 321 AEOrbit%np(icount)=ig 322 AEOrbit%l(icount)=4 323 AEOrbit%kappa(icount)=-5 324 ENDDO 325 ENDIF 326 AEOrbit%norbit=icount 327 AEOrbit%nps=nps 328 AEOrbit%npp=npp 329 AEOrbit%npd=npd 330 AEOrbit%npf=npf 331 AEOrbit%npg=npg 332 IF (AEOrbit%norbit/=input_dataset%norbit) STOP 'Inconsistent number of orbitals!' 333 WRITE(STD_OUT,*) AEOrbit%norbit, ' orbitals will be calculated' 334 335 WRITE(STD_OUT,*)' Below are listed the default occupations ' 336 WRITE(STD_OUT,"(' n l kappa occupancy')") 337 DO io=1,AEOrbit%norbit 338 WRITE(STD_OUT,'(i2,1x,i2,3x,i2,4x,1p,1e15.7)') & 339 & AEOrbit%np(io),AEOrbit%l(io),AEOrbit%kappa(io),AEOrbit%occ(io) 340 ENDDO 341 342 !Corrected occupations (from input dataset) 343 DO io=1,input_dataset%norbit_mod 344 l=input_dataset%orbit_mod_l(io) 345 ip=input_dataset%orbit_mod_n(io) 346 kappa=input_dataset%orbit_mod_k(io) 347 xocc=input_dataset%orbit_mod_occ(io) 348 nfix=nl(ip,kappa) 349 IF (nfix<=0.OR.nfix>AEOrbit%norbit) THEN 350 WRITE(STD_OUT,*) 'error in occupations -- ip,l,kappa,xocc: ',ip,l,kappa,xocc,nfix,AEOrbit%norbit 351 STOP 352 ENDIF 353 AEOrbit%occ(nfix)=xocc 354 END DO 355 356 WRITE(STD_OUT,*) ' Corrected occupations are: ' 357 WRITE(STD_OUT,"(' n l kappa occupancy')") 358 electrons=0.d0 359 DO io=1,AEOrbit%norbit 360 WRITE(STD_OUT,'(i2,1x,i2,3x,i2,4x,1p,1e15.7)') & 361 & AEOrbit%np(io),AEOrbit%l(io),AEOrbit%kappa(io),AEOrbit%occ(io) 362 electrons=electrons+AEOrbit%occ(io) 363 ENDDO 364 ENDIF ! diracrelativistic 365 366 AEPot%q=electrons 367 qf=AEPot%nz-electrons 368 WRITE(STD_OUT,*) 369 WRITE(STD_OUT,*) 'nuclear charge = ', AEPot%nz 370 WRITE(STD_OUT,*) 'electronic charge = ', electrons 371 WRITE(STD_OUT,*) 'net charge = ', qf 372 373 CALL InitSCF(AESCF) 374 IF (scalarrelativistic) CALL Allocate_Scalar_Relativistic(Grid) 375 IF (diracrelativistic) CALL Allocate_Dirac_Relativistic(Grid) 376 377 write(std_out,*) 'Finish SCFatom_Init' ; call flush_unit(std_out) 378 379 DEALLOCATE(nl) 380 381 END SUBROUTINE SCFatom_Init 382 383!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 384!! SCFatom !! 385!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 386 387 SUBROUTINE SCFatom(scftype,lotsofoutput,skip_reading) 388 IMPLICIT NONE 389 CHARACTER(2),INTENT(IN)::scftype 390 LOGICAL,INTENT(IN) :: lotsofoutput 391 LOGICAL,INTENT(IN),OPTIONAL :: skip_reading 392 TYPE(SCFInfo) :: SCFPP 393 394 ! program to perform self-consistency loop 395 396 IF(scftype=='AE') THEN 397 !------------AE calculation-------------- 398 frozencorecalculation=.FALSE. 399 setupfrozencore=.FALSE. 400 PotPtr=>AEPot 401 OrbitPtr=>AEOrbit 402 SCFPtr=>AESCF 403 CALL Orbit_Init(OrbitPtr,PotPtr) 404 CALL Potential_Init(OrbitPtr,PotPtr) 405 406 ELSEIF(scftype=='NC') THEN 407 !------------New Configuration Calculation--------- 408 frozencorecalculation=.FALSE. 409 setupfrozencore=.FALSE. 410 PotPtr=>AEPot 411 OrbitPtr=>AEOrbit 412 SCFPtr=>AESCF 413 IF (PRESENT(skip_reading)) THEN 414 CALL NC_Init(OrbitPtr,PotPtr,skip_reading=skip_reading) 415 ELSE 416 CALL NC_Init(OrbitPtr,PotPtr) 417 ENDIF 418 CALL Potential_Init(OrbitPtr,PotPtr) 419 420 ELSEIF(scftype=='SC') THEN 421 !------------Set Core and Valence in current config.--------- 422 frozencorecalculation=.TRUE. 423 setupfrozencore=.TRUE. 424 CALL CopyPot(AEPot,FCPot) 425 CALL CopyOrbit(AEOrbit,FCOrbit) 426 CALL CopySCF(AESCF,FCSCF) 427 PotPtr=>FCPot 428 OrbitPtr=>FCOrbit 429 SCFPtr=>FCSCF 430 CALL SC_Init(OrbitPtr,PotPtr,SCFPtr) 431 432 ELSEIF(scftype=='FC') THEN 433 !------------Frozen Core Calculation--------- 434 frozencorecalculation=.TRUE. 435 setupfrozencore=.FALSE. 436 PotPtr=>FCPot 437 OrbitPtr=>FCOrbit 438 SCFPtr=>FCSCF 439 CALL FC_Init(OrbitPtr,PotPtr) 440 CALL Potential_Init(OrbitPtr,PotPtr) 441 442 ENDIF 443 444 IF (TRIM(OrbitPtr%exctype)=='EXX'.OR. & 445& TRIM(OrbitPtr%exctype)=='EXXKLI'.OR. & 446 TRIM(OrbitPtr%exctype)=='EXXCS') THEN 447 CALL EXX_SCF(scftype,lotsofoutput,Grid,OrbitPtr,PotPtr,FC,SCFPtr) 448 449 ELSE IF (TRIM(OrbitPtr%exctype)=='EXXOCC') THEN 450 CALL EXXOCC_SCF(scftype,lotsofoutput,Grid,OrbitPtr,PotPtr,FC,SCFPtr) 451 452 ELSE IF (TRIM(OrbitPtr%exctype)=='HF'.or.TRIM(OrbitPtr%exctype)=='HFV') THEN 453 write(std_out,*) 'Just before HF_SCF'; call flush_unit(std_out) 454 CALL HF_SCF(scftype,lotsofoutput,Grid,OrbitPtr,PotPtr,FC,SCFPtr) 455 write(std_out,*) 'Just after HF_SCF'; call flush_unit(std_out) 456 457 ELSE 458 !LDA or GGA 459 CALL LDAGGA_SCF(scftype,lotsofoutput,Grid,OrbitPtr,PotPtr,FC,SCFPtr) 460 ENDIF 461 462 If (HFpostprocess) then 463 write(std_out,*) "******PostProcessing HF*********" 464 CALL InitSCF(SCFPP) 465 call hf_energy_only(Grid,OrbitPtr,PotPtr,SCFPP) 466 endif 467 END SUBROUTINE SCFatom 468 469!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 470 !! Orbit_Init 471 !! From nuclear charge -- generate hydrogenic-like initial wfns 472 !! and densities -- 473 !! fill AEOrbit%wfn, AEOrbit%eig, and AEOrbit%den and AEOrbit%q 474 !! also AEOrbit%otau and AEOrbit%tau 475 !! Note that both den and tau need to be divided by 4 \pi r^2 476 !! Note that tau is only correct for non-relativistic case 477!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 478 SUBROUTINE Orbit_Init(Orbit,Pot) 479 IMPLICIT NONE 480 TYPE(PotentialInfo),INTENT(INOUT) :: Pot 481 TYPE(OrbitInfo),INTENT(INOUT) :: Orbit 482 INTEGER :: i,io,ir,ip,l,nfix,np,kappa 483 REAL(8) :: en0,qcal,qf,rescale,z,small,zeff,xocc,fac 484 REAL(8),allocatable :: dpdr(:),pbr(:) 485 INTEGER :: initialconfig=0 486 487 IF (initialconfig/=0) STOP 'Error in aeatom -- Orbit_Init already called' 488 489 ! calculate initial charge density from hydrogen-like functions 490 ! also initial energies 491 zeff=Pot%nz 492 If (.not.diracrelativistic) then 493 DO io=1,Orbit%norbit 494 np=Orbit%np(io) 495 l=Orbit%l(io) 496 xocc=Orbit%occ(io) 497 Orbit%eig(io)=-(zeff/(np))**2 498 WRITE(STD_OUT,*) io,np,l,xocc,Orbit%eig(io) 499 DO ir=1,Grid%n 500 Orbit%wfn(ir,io)=hwfn(zeff,np,l,Grid%r(ir)) 501 IF (ABS(Orbit%wfn(ir,io))<machine_zero) Orbit%wfn(ir,io)=0.d0 502 ENDDO 503 zeff=zeff-0.5d0*xocc 504 zeff=MAX(zeff,1.d0) 505 ENDDO 506 endif 507 If (diracrelativistic) then 508 DO io=1,Orbit%norbit 509 np=Orbit%np(io) 510 l=Orbit%l(io) 511 kappa=Orbit%kappa(io) 512 xocc=Orbit%occ(io) 513 Orbit%eig(io)=-(zeff/(np))**2 514 WRITE(STD_OUT,*) io,np,l,xocc,Orbit%eig(io) 515 DO ir=1,Grid%n 516 call dirachwfn(np,kappa,zeff,Grid%r(ir),Orbit%eig(io) & 517& ,Orbit%wfn(ir,io),Orbit%lwfn(ir,io)) 518 IF (ABS(Orbit%wfn(ir,io))<machine_zero) Orbit%wfn(ir,io)=0.d0 519 IF (ABS(Orbit%lwfn(ir,io))<machine_zero) Orbit%lwfn(ir,io)=0.d0 520 ENDDO 521 zeff=zeff-0.5d0*xocc 522 zeff=MAX(zeff,1.d0) 523 ENDDO 524 endif 525 526 527 528 ! check charge and rescale 529 allocate(dpdr(Grid%n),pbr(Grid%n)) 530 Orbit%den=0.d0;Orbit%tau=0.d0 531 DO io=1,Orbit%norbit 532 dpdr=0.d0;pbr=0.d0 533 pbr(2:Grid%n)=Orbit%wfn(2:Grid%n,io)/Grid%r(2:Grid%n) 534 CALL derivative(Grid,pbr,dpdr,2,Grid%n) 535 CALL extrapolate(Grid,pbr) 536 CALL extrapolate(Grid,dpdr) 537 l=Orbit%l(io) 538 fac=l*(l+1) 539 Orbit%otau(:,io)=(Grid%r(:)*dpdr(:))**2+fac*(pbr(:))**2 540 xocc=Orbit%occ(io) 541 DO ir=1,Grid%n 542 Orbit%den(ir)=Orbit%den(ir)+xocc*(Orbit%wfn(ir,io)**2) 543 Orbit%tau(ir)=Orbit%tau(ir)+xocc*Orbit%otau(ir,io) 544 If (diracrelativistic) Orbit%den(ir)=Orbit%den(ir) + & 545& xocc*((Orbit%lwfn(ir,io))**2) 546 ENDDO 547 ENDDO 548! Note that kinetic energy density (tau) is in Rydberg units 549! Note that kinetic energy density is only correct for non-relativistic 550! formulation 551 qcal=integrator(Grid,Orbit%den) 552 qf=qcal 553 !WRITE(STD_OUT,*) 'qcal electrons = ',qcal, electrons 554 !rescale density 555 rescale=electrons/qcal 556 Orbit%den(1:Grid%n)=Orbit%den(1:Grid%n)*rescale 557 Orbit%tau(1:Grid%n)=Orbit%tau(1:Grid%n)*rescale 558 559 deallocate(dpdr,pbr) 560 initialconfig=1 561 write(std_out,*) 'completed Orbit_Init '; call flush_unit(std_out) 562 563 END SUBROUTINE Orbit_Init 564 565!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 566 SUBROUTINE NC_Init(Orbit,Pot,skip_reading) 567 TYPE(PotentialInfo),INTENT(INOUT) :: Pot 568 TYPE(OrbitInfo),INTENT(INOUT) :: Orbit 569 LOGICAL,INTENT(IN),OPTIONAL :: skip_reading 570 INTEGER :: i,io,jo,ir,ip,l,nfix,j,norbit_mod,np(5) 571 LOGICAL :: read_new_occ 572 REAL(8) :: en0,qcal,qf,rescale,z,small,zeff,xocc,fac 573 REAL(8),allocatable :: dpdr(:),pbr(:) 574 INTEGER, ALLOCATABLE :: orbit_mod_l(:),orbit_mod_n(:),orbit_mod_k(:) 575 REAL(8), ALLOCATABLE :: orbit_mod_occ(:) 576 577 read_new_occ=.true.;if (present(skip_reading)) read_new_occ=.not.skip_reading 578 579 !----------New Configuration--AE calculation-------------- 580 ! readin revision and normalize the density 581 582 IF (read_new_occ) THEN 583 584 !Read modified occupations from standard input 585 np(1)=Orbit%nps;np(2)=Orbit%npp;np(3)=Orbit%npd;np(4)=Orbit%npf;np(5)=Orbit%npg 586 CALL input_dataset_read_occ(norbit_mod,orbit_mod_l,orbit_mod_n,orbit_mod_k,& 587 & orbit_mod_occ,np,diracrelativistic) 588 589 DO jo=1,norbit_mod 590 nfix=-100 591 DO io=1,Orbit%norbit 592 IF (orbit_mod_n(jo)==Orbit%np(io).AND.orbit_mod_l(jo)==Orbit%l(io).AND.& 593 & ((.NOT.diracrelativistic).OR.orbit_mod_k(jo)==Orbit%kappa(io))) THEN 594 nfix=io 595 EXIT 596 ENDIF 597 ENDDO 598 IF (nfix.LE.0.OR.nfix.GT.Orbit%norbit) THEN 599 WRITE(STD_OUT,*) 'error in occupations -- ip,l,xocc',& 600 orbit_mod_n(jo),orbit_mod_l(jo),orbit_mod_occ(jo),nfix,Orbit%norbit 601 STOP 602 ENDIF 603 Orbit%occ(nfix)=orbit_mod_occ(jo) 604 ENDDO 605 DEALLOCATE(orbit_mod_l) 606 DEALLOCATE(orbit_mod_n) 607 DEALLOCATE(orbit_mod_k) 608 DEALLOCATE(orbit_mod_occ) 609 ENDIF 610 611 electrons=0.d0 612 WRITE(STD_OUT,*) ' Corrected occupations are: ' 613 IF (.NOT.diracrelativistic) THEN 614 WRITE(STD_OUT,"(' n l occupancy')") 615 DO io=1,Orbit%norbit 616 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,1e15.7)') Orbit%np(io),Orbit%l(io),Orbit%occ(io) 617 electrons=electrons+Orbit%occ(io) 618 ENDDO 619 ELSE 620 WRITE(STD_OUT,"(' n l kap occupancy')") 621 DO io=1,Orbit%norbit 622 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,1e15.7)') Orbit%np(io),Orbit%l(io),Orbit%kappa(io),Orbit%occ(io) 623 electrons=electrons+Orbit%occ(io) 624 ENDDO 625 ENDIF 626 Pot%q=electrons 627 qf=Pot%nz-electrons 628 WRITE(STD_OUT,*) 629 WRITE(STD_OUT,*) 'nuclear charge = ' , Pot%nz 630 WRITE(STD_OUT,*) 'electronic charge = ', electrons 631 WRITE(STD_OUT,*) 'net charge = ', qf 632 ! 633 ! 634 ! calculate initial charge density from stored wavefunctions 635 ! also initial energies 636 ! 637 allocate(dpdr(Grid%n),pbr(Grid%n)) 638 Orbit%den=0.d0;Orbit%tau=0.d0 639 DO io=1,Orbit%norbit 640 dpdr=0.d0;pbr=0.d0 641 pbr(2:Grid%n)=Orbit%wfn(2:Grid%n,io)/Grid%r(2:Grid%n) 642 CALL derivative(Grid,pbr,dpdr,2,Grid%n) 643 CALL extrapolate(Grid,pbr) 644 CALL extrapolate(Grid,dpdr) 645 l=Orbit%l(io) 646 fac=l*(l+1) 647 Orbit%otau(:,io)=(Grid%r(:)*dpdr(:))**2+fac*(pbr(:))**2 648 xocc=Orbit%occ(io) 649 DO ir=1,Grid%n 650 Orbit%den(ir)=Orbit%den(ir)+xocc*(Orbit%wfn(ir,io)**2) 651 Orbit%tau(ir)=Orbit%tau(ir)+xocc*Orbit%otau(ir,io) 652 If (diracrelativistic) Orbit%den(ir)=Orbit%den(ir) + & 653& xocc*((Orbit%lwfn(ir,io))**2) 654 ENDDO 655 ENDDO 656! Note that kinetic energy density (tau) is in Rydberg units 657! Note that kinetic energy density is only correct for non-relativistic 658! formulation 659 ! 660 ! check charge 661 ! 662 qcal=integrator(Grid,Orbit%den) 663 qf=qcal 664 !WRITE(STD_OUT,*) 'qcal electrons = ',qcal, electrons 665 ! rescale density 666 rescale=electrons/qcal 667 Orbit%den=Orbit%den*rescale 668 Orbit%tau=Orbit%tau*rescale 669 670 deallocate(dpdr,pbr) 671 END SUBROUTINE NC_Init 672 673!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 674!! SC_Init -- set core states in current configuration 675!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 676 SUBROUTINE SC_Init(Orbit,Pot,SCF) 677 TYPE(PotentialInfo),INTENT(INOUT) :: Pot 678 TYPE(OrbitInfo),INTENT(INOUT) :: Orbit 679 TYPE (SCFInfo),INTENT(INOUT) :: SCF 680 INTEGER :: choosevalence=0 681 682 INTEGER :: io 683 LOGICAL :: noalt=.true. 684 If (choosevalence/=0) then 685 write(std_out,*) 'Error in aeatom -- SC_Init already called' 686 stop 687 endif 688 689 CALL Set_Valence() 690 Call Core_Electron_Report(Orbit,FC,std_out) 691 Call Valence_Electron_Report(Orbit,FC,std_out) 692 693 IF (TRIM(Orbit%exctype)=='EXX') THEN 694 Call Get_FCKinCoul(Grid,Pot,Orbit,FC,SCF) 695 CALL Get_FCEnergy_EXX(Grid,Orbit,FC,SCF) 696 CALL Set_Vxref(AEPot) ! probably not a good idea 697 ELSEIF (TRIM(Orbit%exctype)=='EXXKLI') THEN 698 Call Get_FCKinCoul(Grid,Pot,Orbit,FC,SCF,noalt) 699 CALL Get_FCEnergy_EXX(Grid,Orbit,FC,SCF) 700 ELSEIF (TRIM(Orbit%exctype)=='HF'.or.TRIM(Orbit%exctype)=='HFV') THEN 701 write(std_out,*) ' Completed Setcore '; call flush_unit(std_out) 702 CALL Get_FCEnergy_HF(Grid,Pot,Orbit,FC,SCF) 703 CALL Total_FCEnergy_Report(SCF,std_out) 704 write(std_out,*) ' Completed Setcore '; call flush_unit(std_out) 705 ELSE 706 Call Get_FCKinCoul(Grid,Pot,Orbit,FC,SCF) 707 CALL Get_FCEXC(SCF) 708 CALL Total_FCEnergy_Report(SCF,std_out) 709 ENDIF 710 !CALL Total_Energy_Report(SCF,std_out) 711 ! For EXX, exchange contribution is not known yet. 712 713 choosevalence=1 714 715 END SUBROUTINE SC_Init 716 717!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 718 SUBROUTINE FC_Init(Orbit,Pot) 719 TYPE(PotentialInfo),INTENT(INOUT) :: Pot 720 TYPE(OrbitInfo),INTENT(INOUT) :: Orbit 721 INTEGER :: i,io,jo,ir,ip,l,nfix,j,norbit_mod,np(5) 722 REAL(8) :: en0,qcal,qf,rescale,z,small,zeff,xocc 723 INTEGER, ALLOCATABLE :: orbit_mod_l(:),orbit_mod_n(:),orbit_mod_k(:) 724 REAL(8), ALLOCATABLE :: orbit_mod_occ(:) 725 726 WRITE(STD_OUT,*)' Below are listed the current valence occupations ' 727 IF (.NOT.diracrelativistic) THEN 728 WRITE(STD_OUT,"(' n l occupancy')") 729 DO io=1,Orbit%norbit 730 IF (.NOT.Orbit%iscore(io)) WRITE(STD_OUT,'(i2,1x,i2,4x,1p,1e15.7)') & 731& Orbit%np(io),Orbit%l(io),Orbit%occ(io) 732 ENDDO 733 ELSE 734 WRITE(STD_OUT,"(' n l kap occupancy')") 735 DO io=1,Orbit%norbit 736 IF (.NOT.Orbit%iscore(io)) WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,1e15.7)') & 737& Orbit%np(io),Orbit%l(io),Orbit%kappa(io),Orbit%occ(io) 738 ENDDO 739 ENDIF 740 741 !Read modified occupations from standard input 742 np(1)=Orbit%nps;np(2)=Orbit%npp;np(3)=Orbit%npd;np(4)=Orbit%npf;np(5)=Orbit%npg 743 CALL input_dataset_read_occ(norbit_mod,orbit_mod_l,orbit_mod_n,orbit_mod_k,& 744& orbit_mod_occ,np,diracrelativistic) 745 746 DO jo=1,norbit_mod 747 nfix=-100 748 DO io=1,Orbit%norbit 749 IF (orbit_mod_n(jo)==Orbit%np(io).AND.orbit_mod_l(jo)==Orbit%l(io).AND.& 750& ((.NOT.diracrelativistic).OR.orbit_mod_k(jo)==Orbit%kappa(io))) THEN 751 nfix=io 752 EXIT 753 ENDIF 754 ENDDO 755 IF (nfix.LE.0.OR.nfix.GT.Orbit%norbit) THEN 756 WRITE(STD_OUT,*) 'error in occupations -- ip,l,xocc',& 757 orbit_mod_n(jo),orbit_mod_l(jo),orbit_mod_occ(jo),nfix,Orbit%norbit 758 STOP 759 ENDIF 760 Orbit%occ(nfix)=orbit_mod_occ(jo) 761 ENDDO 762 DEALLOCATE(orbit_mod_l) 763 DEALLOCATE(orbit_mod_n) 764 DEALLOCATE(orbit_mod_k) 765 DEALLOCATE(orbit_mod_occ) 766 767 FC%zvale=0.d0 768 WRITE(STD_OUT,*) ' Corrected occupations are: ' 769 IF (.NOT.diracrelativistic) THEN 770 WRITE(STD_OUT,"(' n l occupancy')") 771 DO io=1,Orbit%norbit 772 If (.not.Orbit%iscore(io))then 773 WRITE(STD_OUT,'(i2,1x,i2,4x,1p,1e15.7)') & 774& Orbit%np(io),Orbit%l(io),Orbit%occ(io) 775 FC%zvale=FC%zvale+Orbit%occ(io) 776 Endif 777 ENDDO 778 ELSE 779 WRITE(STD_OUT,"(' n l kap occupancy')") 780 DO io=1,Orbit%norbit 781 If (.not.Orbit%iscore(io))then 782 WRITE(STD_OUT,'(i2,1x,i2,2x,i2,4x,1p,1e15.7)') & 783& Orbit%np(io),Orbit%l(io),Orbit%kappa(io),Orbit%occ(io) 784 FC%zvale=FC%zvale+Orbit%occ(io) 785 Endif 786 ENDDO 787 ENDIF 788 789 electrons=FC%zvale+FC%zcore 790 791 ! 792 ! calculate initial charge density from stored wavefunctions 793 ! also initial energies 794 ! 795 FC%valeden=0; FC%valetau=0 796 DO io=1,Orbit%norbit 797 IF (.NOT.Orbit%iscore(io)) THEN 798 xocc=Orbit%occ(io) 799 DO ir=1,Grid%n 800 FC%valeden(ir)=FC%valeden(ir)+xocc*(Orbit%wfn(ir,io)**2) 801 FC%valetau(ir)=FC%valetau(ir)+xocc*Orbit%otau(ir,io) 802 If (diracrelativistic) FC%valeden(ir)=FC%valeden(ir) + & 803& xocc*((Orbit%lwfn(ir,io))**2) 804 805 ENDDO 806 ENDIF 807 ENDDO 808 ! 809 ! check charge ! rescale density 810 811 ! 812 qcal=integrator(Grid,FC%valeden) 813 if (ABS(qcal)<small0.OR.electrons<small0) then 814 FC%valeden=0 815 else 816 rescale=FC%zvale/qcal 817 FC%valeden=FC%valeden*rescale 818 FC%valetau=FC%valetau*rescale 819 endif 820 821 Orbit%den=FC%valeden+FC%coreden 822 WRITE(STD_OUT,*) 823 WRITE(STD_OUT,*) 'core charge = ' , FC%zcore 824 WRITE(STD_OUT,*) 'valence charge = ', FC%zvale 825 WRITE(STD_OUT,*) 'net charge = ', Pot%nz-FC%zcore-FC%zvale 826 827 END SUBROUTINE FC_Init 828 829!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 830 !! Potential_Init 831 !! Generate Potptr%rv, rvh, rvx given AEOrbit%wfn, AEOrbit%den, AEOrbit%q 832!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 833 SUBROUTINE Potential_Init(Orbit,Pot) 834 IMPLICIT NONE 835 TYPE(OrbitInfo),INTENT(IN) :: Orbit 836 TYPE(PotentialInfo),INTENT(INOUT) :: Pot 837 INTEGER :: i,io,ir,xocc,ip,l,nfix,j,np 838 REAL(8) :: en0,qcal,qf,rescale,z,small,zeff,ecoul,v0,etxc,eex 839 LOGICAL :: success 840 841 842 CALL poisson(Grid,Pot%q,Orbit%den,Pot%rvh,ecoul,v00=v0) 843 write(std_out,*) 'In Potential_Init', Pot%q,ecoul; call flush_unit(std_out) 844 845 END SUBROUTINE Potential_Init 846 847 848!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 849 ! Set Valence Orbitals 850 ! 851!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 852 SUBROUTINE Set_Valence() 853 854 INTEGER :: io,n,ok 855 856 call InitFC(FC,Grid%n) 857 858 FC%zcore=0; FC%zvale=0 859 FC%coreden=0; FC%valeden=0 860 FC%coretau=0; FC%valetau=0 861 862 DO io=1,FCOrbit%norbit 863 WRITE(STD_OUT,'(3i5,1p,2e15.7)') io,FCOrbit%np(io),FCOrbit%l(io),& 864& FCOrbit%occ(io),FCOrbit%eig(io) 865 866 FCOrbit%iscore(io)=input_dataset%orbit_iscore(io) 867 868 IF (FCOrbit%iscore(io)) THEN 869 FC%zcore=FC%zcore+FCOrbit%occ(io) 870 FC%coreden=FC%coreden+FCOrbit%occ(io)*(FCOrbit%wfn(:,io))**2 871 FC%coretau=FC%coretau+FCOrbit%occ(io)*FCOrbit%otau(:,io) 872 If (diracrelativistic) FC%coreden=FC%coreden + & 873& FCOrbit%occ(io)*((FCOrbit%lwfn(:,io))**2) 874 875 ENDIF 876 IF (.NOT.FCOrbit%iscore(io)) THEN 877 878 FC%zvale=FC%zvale+FCOrbit%occ(io) 879 FC%valeden=FC%valeden+FCOrbit%occ(io)*(FCOrbit%wfn(:,io))**2 880 FC%valetau=FC%valetau+FCOrbit%occ(io)*FCOrbit%otau(:,io) 881 If (diracrelativistic) FC%valeden=FC%valeden + & 882& FCOrbit%occ(io)*((FCOrbit%lwfn(:,io))**2) 883 ENDIF 884 ENDDO 885 886 write(std_out,*) 'Returning from SetValence' 887 write(std_out,*) 'Core electrons', FC%zcore,integrator(Grid,FC%coreden) 888 write(std_out,*) 'Vale electrons', FC%zvale,integrator(Grid,FC%valeden) 889 890 write(std_out,*) 'Coretau:',FC%coretau(1:20) 891 write(std_out,*) 'Valetau:',FC%valetau(1:20) 892 END SUBROUTINE Set_Valence 893 894! dump and load subroutines need to be updated!!!! 895 SUBROUTINE dump_aeatom(Fn,Grid,AEOrbit,AEPot,AESCF,FCOrbit,FCPot,FCSCF,FC) 896 Character(*), INTENT(IN) :: Fn 897 Type(GridInfo), INTENT(IN) :: Grid 898 Type(OrbitInfo), INTENT(IN) :: AEOrbit,FCOrbit 899 Type(PotentialInfo), INTENT(IN) :: AEPot,FCPot 900 Type(SCFInfo), INTENT(IN) :: AESCF,FCSCF 901 Type(FCInfo), INTENT(IN) :: FC 902 903 INTEGER, parameter :: ifo=15 904 INTEGER :: i,io,n 905 906 Write(STD_OUT,*) 'Note that dump has been called; code not completely checked' 907 908 open(ifo,file=Fn,form='formatted',status='replace') 909 write(std_out,*) 'Creating or replacing dump file', Fn 910 911 write(ifo,*) 'ATOMPAW' ! keyword 912 913 ! grid info 914 write(ifo,*) Grid%TYPE,Grid%n,Grid%h 915 n=Grid%n 916 write(ifo,*)(Grid%r(i),i=1,n) 917 if (usingloggrid(Grid)) then 918 write(ifo,*)(Grid%drdu(i),Grid%pref(i),Grid%rr02(i),i=1,n) 919 endif 920 921 ! atomdata fixed constants 922 write(ifo,*) frozencorecalculation,setupfrozencore,scalarrelativistic,& 923& finitenucleus,gaussianshapefunction,besselshapefunction,& 924& ColleSalvetti 925 926 ! Orbit info 927 write(ifo,*) AEOrbit%exctype,AEOrbit%nps, AEOrbit%npp, AEOrbit%npd ,& 928& AEOrbit%npf, AEOrbit%npg, AEOrbit%norbit 929 write(ifo,*) (AEOrbit%np(io),AEOrbit%l(io),AEOrbit%iscore(io),& 930& AEOrbit%eig(io),AEOrbit%occ(io),io=1,AEOrbit%norbit) 931 write(ifo,*) ((AEOrbit%wfn(i,io),i=1,n),io=1,AEOrbit%norbit) 932 write(ifo,*) (AEOrbit%den(i),i=1,n) 933 write(ifo,*) (AEOrbit%tau(i),i=1,n) 934 935 write(ifo,*) FCOrbit%exctype,FCOrbit%nps, FCOrbit%npp, FCOrbit%npd ,& 936& FCOrbit%npf, FCOrbit%npg, FCOrbit%norbit 937 write(ifo,*) (FCOrbit%np(io),FCOrbit%l(io),FCOrbit%iscore(io),& 938& FCOrbit%eig(io),FCOrbit%occ(io),io=1,FCOrbit%norbit) 939 write(ifo,*) ((FCOrbit%wfn(i,io),i=1,n),io=1,FCOrbit%norbit) 940 write(ifo,*) (FCOrbit%den(i),i=1,n) 941 942 ! Pot info 943 write(ifo,*) AEPot%nz,AEPot%sym,AEPot%q,AEPot%v0,AEPot%v0p 944 write(ifo,*) (AEPot%rv(i),AEPot%rvn(i),AEPot%rvh(i),AEPot%rvx(i),i=1,n) 945 946 write(ifo,*) FCPot%nz,FCPot%sym,FCPot%q,FCPot%v0,FCPot%v0p 947 write(ifo,*) (FCPot%rv(i),FCPot%rvn(i),FCPot%rvh(i),FCPot%rvx(i),i=1,n) 948 949 ! SCFinfo 950 write(ifo,*) AESCF%iter,AESCF%delta,AESCF%eone,AESCF%ekin,& 951& AESCF%estatic,AESCF%ecoul,AESCF%eexc,AESCF%oepcs,AESCF%etot,& 952& AESCF%valekin,AESCF%valecoul,AESCF%valeexc,AESCF%corekin,AESCF%evale 953 954 write(ifo,*) FCSCF%iter,FCSCF%delta,FCSCF%eone,FCSCF%ekin,& 955& FCSCF%estatic,FCSCF%ecoul,FCSCF%eexc,FCSCF%oepcs,FCSCF%etot,& 956& FCSCF%valekin,FCSCF%valecoul,FCSCF%valeexc,FCSCF%corekin,FCSCF%evale 957 958 ! FCinfo 959 write(ifo,*) FC%zvale,FC%zcore 960 write(ifo,*) (FC%coreden(i),FC%valeden(i),i=1,n) 961 write(ifo,*) (FC%coretau(i),FC%valetau(i),i=1,n) 962 963 If (TRIM(AEOrbit%exctype)=='EXX'.or.TRIM(AEOrbit%exctype)=='EXXKLI') & 964& Call EXXdump(Grid,AEOrbit,ifo) 965 If (TRIM(AEOrbit%exctype)=='HF'.or.TRIM(AEOrbit%exctype)=='HFV') & 966& Call HFdump(Grid,AEOrbit,ifo) 967 968 close(ifo) 969 970 write(std_out,*) 'Closing dump file' 971 END SUBROUTINE dump_aeatom 972 973 SUBROUTINE load_aeatom(Fn,Grid,AEOrbit,AEPot,AESCF,FCOrbit,FCPot,FCSCF,FC) 974 Character(*), INTENT(IN) :: Fn 975 Type(GridInfo), INTENT(OUT) :: Grid 976 Type(OrbitInfo), INTENT(OUT) :: AEOrbit,FCOrbit 977 Type(PotentialInfo), INTENT(OUT) :: AEPot,FCPot 978 Type(SCFInfo), INTENT(OUT) :: AESCF,FCSCF 979 Type(FCInfo), INTENT(OUT) :: FC 980 981 INTEGER, parameter :: ifo=15 982 INTEGER :: n,norbit,nps,npp,npd,npf,npg,i,io 983 CHARACTER(132) :: inputfile 984 985 Write(STD_OUT,*) 'Note that load has been called; code not completely checked' 986 987 open(ifo,file=Fn,form='formatted',status='old') 988 write(std_out,*) 'Reading dump file ', Fn 989 990 read(ifo,*) inputfile(1:7) 991 if (inputfile(1:7)=='ATOMPAW') then 992 write(std_out,*) 'File seems to be fine' 993 else 994 write(std_out,*) 'File is not correct -- program will stop' 995 stop 996 endif 997 998 ! grid info 999 read(ifo,*) Grid%TYPE,Grid%n,Grid%h 1000 n=Grid%n 1001 Allocate(Grid%r(n)) 1002 read(ifo,*)(Grid%r(i),i=1,n) 1003 if (usingloggrid(Grid)) then 1004 Allocate(Grid%drdu(n),Grid%pref(n),Grid%rr02(n)) 1005 read(ifo,*)(Grid%drdu(i),Grid%pref(i),Grid%rr02(i),i=1,n) 1006 endif 1007 1008 ! atomdata fixed constants 1009 read(ifo,*) frozencorecalculation,setupfrozencore,scalarrelativistic ,& 1010& finitenucleus,gaussianshapefunction,besselshapefunction,ColleSalvetti 1011 1012 ! Orbit info 1013 read(ifo,*) exctype,nps,npp,npd,npf,npg,norbit 1014 Call InitOrbit(AEOrbit,norbit,n,exctype) 1015 AEOrbit%nps=nps;AEOrbit%npp=npp;AEOrbit%npd=npd;AEOrbit%npf=npf;AEOrbit%npg=npg 1016 read(ifo,*) (AEOrbit%np(io),AEOrbit%l(io),AEOrbit%iscore(io),& 1017& AEOrbit%eig(io),AEOrbit%occ(io),io=1,AEOrbit%norbit) 1018 read(ifo,*) ((AEOrbit%wfn(i,io),i=1,n),io=1,AEOrbit%norbit) 1019 read(ifo,*) (AEOrbit%den(i),i=1,n) 1020 read(ifo,*) (AEOrbit%tau(i),i=1,n) 1021 1022 read(ifo,*) exctype,nps,npp,npd,npf,npg,norbit 1023 Call InitOrbit(FCOrbit,norbit,n,exctype) 1024 FCOrbit%nps=nps;FCOrbit%npp=npp;FCOrbit%npd=npd;FCOrbit%npf=npf;FCOrbit%npg=npg 1025 read(ifo,*) (FCOrbit%np(io),FCOrbit%l(io),FCOrbit%iscore(io),& 1026& FCOrbit%eig(io),FCOrbit%occ(io),io=1,FCOrbit%norbit) 1027 read(ifo,*) ((FCOrbit%wfn(i,io),i=1,n),io=1,FCOrbit%norbit) 1028 read(ifo,*) (FCOrbit%den(i),i=1,n) 1029 1030 exctype=AEOrbit%exctype 1031 IF(TRIM(exctype)/='EXX'.and.TRIM(exctype)/='EXXCS'.and.TRIM(exctype)/='EXXKLI') CALL initexch 1032 1033 ! Pot info 1034 Call InitPot(AEPot,n) 1035 read(ifo,*) AEPot%nz,AEPot%sym,AEPot%q,AEPot%v0,AEPot%v0p 1036 read(ifo,*) (AEPot%rv(i),AEPot%rvn(i),AEPot%rvh(i),AEPot%rvx(i),i=1,n) 1037 1038 Call InitPot(FCPot,n) 1039 read(ifo,*) FCPot%nz,FCPot%sym,FCPot%q,FCPot%v0,FCPot%v0p 1040 read(ifo,*) (FCPot%rv(i),FCPot%rvn(i),FCPot%rvh(i),FCPot%rvx(i),i=1,n) 1041 1042 ! SCFinfo 1043 Call InitSCF(AESCF) 1044 read(ifo,*) AESCF%iter,AESCF%delta,AESCF%eone,AESCF%ekin,& 1045& AESCF%estatic,AESCF%ecoul,AESCF%eexc,AESCF%oepcs,AESCF%etot,& 1046& AESCF%valekin,AESCF%valecoul,AESCF%valeexc,AESCF%corekin,AESCF%evale 1047 1048 Call InitSCF(FCSCF) 1049 read(ifo,*) FCSCF%iter,FCSCF%delta,FCSCF%eone,FCSCF%ekin,& 1050& FCSCF%estatic,FCSCF%ecoul,FCSCF%eexc,FCSCF%oepcs,FCSCF%etot,& 1051& FCSCF%valekin,FCSCF%valecoul,FCSCF%valeexc,FCSCF%corekin,FCSCF%evale 1052 1053 ! FCinfo 1054 call InitFC(FC,n) 1055 read(ifo,*) FC%zvale,FC%zcore 1056 read(ifo,*) (FC%coreden(i),FC%valeden(i),i=1,n) 1057 read(ifo,*) (FC%coretau(i),FC%valetau(i),i=1,n) 1058 1059 If (TRIM(AEOrbit%exctype)=='EXX'.or.TRIM(AEOrbit%exctype)=='EXXKLI') & 1060& Call EXXload(Grid,AEOrbit,ifo) 1061 If (TRIM(AEOrbit%exctype)=='HF'.or.TRIM(AEOrbit%exctype)=='HFV') & 1062& Call HFload(Grid,AEOrbit,ifo) 1063 1064 close(ifo) 1065 1066 write(std_out,*) 'Closing load file' 1067 write(std_out,*) 'Load completed normally' 1068 1069 END SUBROUTINE load_aeatom 1070 1071END MODULE AEatom 1072