DIMENSION PSTOR(0:1400,11),FOCK(1400,11),PHI(0:1400)
dimension RHO(1400),NOCC(11),EXH(11),AML(11),FACTORIAL(0:36)
dimension PSIIN(1400),PSIOUT(0:1400)
COMMON /GENE1/DR,NR/GENE2/E2/GENE3/NSTATE/WAVE/ZSTAR,ABOHR
//KIKAKU/NEE,PFLAG,NIT,VENTOT,VEETOT,VEXTOT,KTOT,ELAST,ETOT
//MITUDO/MIX/KYOU1/ZE2,HBM/KYOU2/NOCC/KYOU3/AML/HIDOUJI/ISTATE,E
//JIGEN1/PSTOR/JIGEN2/RHO/JIGEN3/PHI/JIGEN4/FOCK/JIGEN5/FACTORIAL
//JIGEN6/PSIIN,PSIOUT/JIGEN7/ID/JIGEN8/EXH
REAL*8 PSTOR,FOCK,PHI,RHO, EXH,HBM,E2,ABOHR,ZNC,ZE2,DR,RMAX,RWS,
/ZSTAR,MIX,E,ENUE,KTOT,VENTOT,VEETOT,VEXTOT,KEN,VEN,VEE,VEX,ANG,
/VCENT,PM,R,PZ,PZ2,XTF,FACTORIAL,PSIIN,PSIOUT,AMU,MASS,DENS,PAI,
/RHO0,LACAL,LACCL,MNAD,OMA,ELAST,ETOT,PERVOL,VOLUNI,VOLPRI,
/UNIATOM,PRIATOM
INTEGER*4 NSTATE,NRMAX,NOCC,NR,PFLAG,AML,NIT,SIMASU,NEE,STATE,
/ISTATE,ITERAT,LAT,SURU1,SURU2
CHARACTER*14 ID(11)
CHARACTER*14 FILENAME(2)
C**Start Data ********************* ** ****************
NRMAX=1400
NSTATE=11
HBM=7.6359
E2=14.409
ABOHR=HBM/E2
AMU=1.66053D-24
PAI=3.141592654
ANG=1.D-8
RHO0=0.
C**** Input file-name ************************************
WRITE(6,*)'Do you save HFS Potentials ? 1.Y 2.N'
READ(5,*)SURU1
IF (SURU1 .EQ. 1) THEN
READ (5,10)FILENAME(1)
ENDIF
WRITE(6,*)'the density for Muffin-tin Potentials ? 1.Y 2.N'
READ(5,*)SURU2
IF (SURU2 .EQ. 1) THEN
READ(5,10)FILENAME(2)
ENDIF
10 FORMAT(A14)
C****Input Ion Perameter ******************
WRITE(6,*)'Input Atomic Number Z'
READ(5,*)ZNC
ZE2=ZNC*E2
20 WRITE(6,30)ID(1)
30 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=2)')
READ(5,*)NOCC(1)
IF((NOCC(1) .LT. 0) .OR. (NOCC(1) .GT. 2)) GOTO 20
40 WRITE(6,50)ID(2)
50 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=2)')
READ(5,*)NOCC(2)
IF((NOCC(2) .LT. 0) .OR. (NOCC(2) .GT. 2)) GOTO 40
60 WRITE(6,70)ID(3)
70 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=6)')
READ(5,*)NOCC(3)
IF((NOCC(3) .LT. 0) .OR. (NOCC(3) .GT. 6)) GOTO 60
80 WRITE(6,90)ID(4)
90 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=2)')
READ(5,*)NOCC(4)
IF((NOCC(4) .LT. 0) .OR. (NOCC(4) .GT. 2)) GOTO 80
100 WRITE(6,110)ID(5)
110 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=6)')
READ(5,*)NOCC(5)
IF((NOCC(5) .LT. 0) .OR. (NOCC(5) .GT. 6)) GOTO 100
120 WRITE(6,130)ID(6)
130 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=10)')
READ(5,*) NOCC(6)
IF ((NOCC(6) .LT. 0) .OR. (NOCC(6) .GT. 10)) GOTO 120
140 WRITE(6,150)ID(7)
150 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=2)')
READ(5,*)NOCC(7)
IF((NOCC(7) .LT. 0) .OR. (NOCC(7) .GT. 2)) GOTO 140
160 WRITE(6,170)ID(8)
170 FORMAT(1H ,A14,'the number of electrons in the orbit (>=0,<=6)')
READ(5,*)NOCC(8)
IF((NOCC(8) .LT. 0) .OR. (NOCC(8) .GT. 6)) GOTO 160
180 WRITE(6,190)ID(9)
190 FORMAT(1H ,A14,'the number of electrons in the orbit(>=0,<=10)')
READ (5,*) NOCC(9)
IF((NOCC(9) .LT. 0) .OR. (NOCC(9) .GT. 10)) GOTO 180
200 WRITE(6,210)ID(10)
210 FORMAT(1H ,A14,'the number of electrons in the orbit(>=0.<=14)')
READ(5,*)NOCC(10)
IF((NOCC(10) .LT. 0) .OR. (NOCC(10) .GT. 14)) GOTO 200
220 WRITE(6,230)ID(11)
230 FORMAT(1H ,A14,'the number of electrons in the orbit(>=0,<=2)')
READ(5,*)NOCC(11)
IF((NOCC(11) .LT. 0) .OR. (NOCC(11) .GT. 2)) GOTO 220
WRITE(6,*)'latice parameter a,c (A)'
READ(5,*) LACAL,LACCL
WRITE(6,*)'Muffin-tin radius'
READ(5,*) MNAD
WRITE(6,*)'the crystal structure'
WRITE(6,*)'1. FCC 2. BCC 3. DIAMOND 4. HCP 5. PARTICULAR'
READ(5,*) LAT
IF (LAT .EQ. 1) THEN
VOLUNI=LACAL**3.
UNIATOM=4.
VOLPRI=VOLUNI/UNIATOM
PRIATOM=1.
PERVOL=VOLPRI/PRIATOM
ELSEIF (LAT .EQ. 2) THEN
VOLUNI=LACAL**3.
UNIATOM=2.
VOLPRI=VOLUNI/UNIATOM
PRIATOM=1.
PERVOL=VOLPRI/PRIATOM
ELSEIF (LAT .EQ. 3) THEN
VOLUNI=LACAL**3.
UNIATOM=8.
VOLPRI=VOLUNI/UNIATOM
PRIATOM=1.
PERVOL=VOLPRI/PRIATOM
ELSEIF (LAT .EQ. 4) THEN
VOLPRI=(LACAL**2.)*DSIN(120.*PAI/180.)*LACCL
PRIATOM=2.
PERVOL=VOLPRI/PRIATOM
ELSE
WRITE(6,*)'the unit sell volume (A^3)'
READ(5,*) VOLPRI
WRITE(6,*)'the number of atoms in the unit sell'
READ(5,*) PRIATOM
PERVOL=VOLPRI/PRIATOM
ENDIF
RWS=((PERVOL*3./4./PAI)**(1./3.))
240 WRITE(6,*)'Wigner-Seitz radius ',RWS,' (A)'
WRITE(6,*)'The radius is equal to about 4(A) or 5(A)'
WRITE(6,*)'That is longer than Wigner-Seitz radius at least'
WRITE(6,*)'The limit is 7(A) for atomic potential'
WRITE(6,*)'pith ( 0.005 A - 0.3 A)'
READ(5,*) DR
WRITE(6,*)'the radius of the shell (A)'
READ(5,*)RMAX
NR=IDNINT(RMAX/DR)
IF(NR .LE. NRMAX)GOTO 270
WRITE(6,250)NR
250 FORMAT(1H ,'the number of the ratice spot ',I5)
WRITE(6,260)NRMAX
260 FORMAT(1H ,'the maxit number of the ratice spot ',I5)
WRITE(6,*)'retry'
GOTO 240
270 WRITE(6,*)'the number of the iteration'
READ(5,*)ITERAT
NEE=0
DO 280 STATE=1,NSTATE
NEE=NEE+NOCC(STATE)
280 CONTINUE
WRITE(6,*)'the summention of electrons ',NEE
WRITE(6,*)'******* HANG ON. I AM WORKING ON IT !! ******'
C************************************************************
ZSTAR=ZNC
C****************************************************
C WRITE(6,*)'CALL NORMALIZATION'
CALL NORMWAVE
MIX=1.
PFLAG=0
C****************************
C**********************
C WRITE(6,*)'CALL ENERGY'
CALL ENERGY
ZSTAR=-ZNC*(VENTOT+VEETOT+VEXTOT)/(2.*KTOT)
C*************************************************
CALL NORMWAVE
WRITE(6,290) ZNC
290 FORMAT(1H ,'atomic number Z= ',F5.1)
WRITE(6,300)NEE
300 FORMAT(1H ,'the number of electorons = ',I5)
WRITE(6,310)DR
310 FORMAT(1H ,'Step = ',F10.5)
WRITE(6,320)DR*DBLE(NR)
320 FORMAT(1H ,'maxit radius = ',F15.5)
WRITE(6,330)ZSTAR
330 FORMAT(1H ,'Zeft = ',F20.8)
PFLAG=-1
NIT=0
ELAST=0.
C***************************
C****************************************
C WRITE(6,*)'CALL ENERGY'
CALL ENERGY
WRITE(6,340)ETOT
340 FORMAT(1H ,'Total energy = ',E15.8)
MIX=0.5
C********Hartree Fock****************************************
DO 360 NIT=1,ITERAT
DO 350 ISTATE=1,NSTATE
IF(NOCC(ISTATE) .EQ. 0)GOTO 350
E=EXH(ISTATE)
IF(E .GT. 0.)E=-10.
C ***************************************************
C WRITE(6,*)'ISTATE,CALL IMHOMO ',ISTATE
CALL IMHOMO
350 CONTINUE
C ***************************
C *****************************************
CALL ENERGY
360 CONTINUE
C***** Save *********************** *****
IF (SURU1 .EQ. 1) THEN
GOTO 370
ELSEIF (SURU2 .EQ. 1) THEN
GOTO 449
ELSE
GOTO 450
ENDIF
370 WRITE(6,*)'****** Now I am saving HFS Potentials *****'
OPEN(15,FILE=FILENAME(1),STATUS='unknown')
C WRITE(15,380)DR,NR
WRITE(15,*) DR,',',NR
CC2022.1.23 380 FORMAT(D15.8,I7)
CCCCCCC 2022.1.23CCCCCCCCCCCC start
C WRITE(15,390)(((-ZE2/(DR*DBLE(I)))+PHI(I)-(3.*((3./(2.**5.)
C / /(PAI**2.)*RHO(I)/((DBLE(I)*DR)**2.))**(1./3.))*E2)),I=1,NR)
C 390 FORMAT(5D15.8)
CCCCCCC 2022.1.23CCCCCCCCCCCC end
do 991 I=1,NR
WRITE(15,*)i,',',-((-ZE2/(DR*DBLE(I)))+PHI(I)-(3.*((3./(2.**5.)
/ /(PAI**2.)*RHO(I)/((DBLE(I)*DR)**2.))**(1./3.))*E2))
991 continue
CLOSE(15)
WRITE(6,*)'*** 0.K ***'
IF (SURU2 .NE. 1) GOTO 450
449 WRITE(6,*)'*** Now I am saving the density for muffin-tin *** '
OPEN(14,FILE=FILENAME(2),STATUS='unknown')
WRITE(14,373) ZNC,MNAD,LAT
373 FORMAT(2D15.8,I7)
WRITE(14,374) PERVOL,LACAL,LACCL
374 FORMAT(3D15.8)
WRITE(14,371) DR,RMAX,NR+1
371 FORMAT(2D15.8,I7)
WRITE(14,372) RHO0,(RHO(I),I=1,NR)
C WRITE(14,372) (PHI(I),I=1,NR)
372 FORMAT(5D15.8)
CLOSE(14)
WRITE(6,*)'*** O.K ***'
450 STOP
END
C***** Subroutin*******************************
C****'set start value' Subroutin ***********************
C **** parameter of each orbit **********************
BLOCK DATA
DIMENSION AML(11)
COMMON /KYOU3/AML/JIGEN7/ID
INTEGER*4 AML
CHARACTER*14 ID(11)
DATA (AML(I),I=1,11)/0,0,1,0,1,2,0,1,2,3,0/,ID/'1s','2s',
/'2p','3s','3p','3d','4s','4p','4d','4f','5s'/
END
C******************************************************************
SUBROUTINE NORMWAVE
dimension NOCC(11),PSTOR(0:1400,11)
COMMON /GENE1/DR,NR/GENE3/NSTATE/WAVE/ZSTAR,ABOHR/KYOU2/NOCC
//JIGEN1/PSTOR
REAL*8 ZSTAR,DR,PSTOR,ABOHR,RSTAR,ESTAR2,NORM,ESTAR3,ESTAR4,
/ESTAR5,DPSTOR,DPLOOG
INTEGER*4 NOCC,NR,NSTATE,STATE
C WRITE(6,*)'Normalizatiion'
DO 110 J=0,NR
RSTAR=DBLE(J)*DR*ZSTAR/ABOHR
IF((RSTAR/2.) .GE. 705.4) THEN
ESTAR2=0.
ELSE
ESTAR2=DEXP(-RSTAR/2.)
ENDIF
IF((RSTAR/3.) .GE. 705.4)THEN
ESTAR3=0.
ELSE
ESTAR3=DEXP(-RSTAR/3.)
ENDIF
IF((RSTAR/4.) .GE. 705.4) THEN
ESTAR4=0.
ELSE
ESTAR4=DEXP(-RSTAR/4.)
ENDIF
IF((RSTAR/5.) .GE. 705.4) THEN
ESTAR5=0.
ELSE
ESTAR5=DEXP(-RSTAR/5.)
ENDIF
C WRITE(6,*)'J,ESTAR2,ESTAR3',J,ESTAR2,ESTAR3
C WRITE(6,*)'J,ESTAR4,ESTAR5',J,ESTAR4,ESTAR5
IF(NOCC(1) .EQ. 0)GOTO 10
c write(*,*) 'AAAAAAAAA'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR2 .EQ. 0.))THEN
C write(*,*) 'lklkasjdoirjwr'
PSTOR(J,1)=0.
ELSE
C write(*,*) 'KOJHWJFHWHRF',Rstar,Estar2
DPSTOR=LOG10(RSTAR)+2.*LOG10(ESTAR2)
C write(*,*) 'KOJHWJFHWHRF'
IF(DPSTOR .LT. ((-307.)+LOG10(0.9)))THEN
PSTOR(J,1)=0.
ELSE
PSTOR(J,1)=RSTAR*(ESTAR2**2.)
C write(*,*) 'RSTAR,ESTAR2',RSTAR,ESTAR2
C write(*,*) 'j,pstor',j,pstor(j,1)
ENDIF
c write(*,*) '1213'
ENDIF
10 IF(NOCC(2) .EQ. 0) GOTO 20
c write(*,*) 'BBBBBBBBBBBBB'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR2 .EQ. 0.) .OR.
/ ((2.-RSTAR) .EQ. 0.)) THEN
PSTOR(J,2)=0.
ELSE
DPSTOR=LOG10(ABS(2.-RSTAR))+LOG10(RSTAR)+LOG10(ESTAR2)
IF(DPSTOR .LT. ((-307.) +LOG10(4.31)))THEN
PSTOR(J,2)=0.
ELSE
PSTOR(J,2)=(2.-RSTAR)*RSTAR*ESTAR2
ENDIF
ENDIF
20 IF(NOCC(3) .EQ. 0) GOTO 30
c write(*,*) 'CCCCCCCCCCCCCCCCC'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR2 .EQ. 0.)) THEN
PSTOR(J,3)=0.
ELSE
DPSTOR=2.*LOG10(RSTAR)+LOG10(ESTAR2)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,3)=0.
ELSE
PSTOR(J,3)=(RSTAR**2.)*ESTAR2
ENDIF
ENDIF
30 IF(NOCC(4) .EQ. 0) GOTO 40
c write(*,*) 'DDDDDDDDDDDDDD'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR3 .EQ. 0.) .OR.
/ ((27.-18.*RSTAR+2.*(RSTAR**2.)) .EQ. 0.))THEN
PSTOR(J,4)=0.
ELSE
DPSTOR=LOG10(ABS(27.-18.*RSTAR+2.*(RSTAR**2.)))
/ +LOG10(RSTAR)+LOG10(ESTAR3)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,4)=0.
ELSE
PSTOR(J,4)=(27.-18.*RSTAR+2.*(RSTAR**2.))*RSTAR*ESTAR3
ENDIF
ENDIF
c write(*,*) 'EEEEEEEEEEEEEE'
40 IF(NOCC(5) .EQ. 0) GOTO 50
c write(*,*) 'FFFFFFFFFFFFFF'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR3 .EQ. 0.).OR.
/ ((6.-RSTAR).EQ. 0.))THEN
PSTOR(J,5)=0.
ELSE
DPSTOR=LOG10(ABS(6.-RSTAR))
/ +2.*LOG10(RSTAR)+LOG10(ESTAR3)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,5)=0.
ELSE
PSTOR(J,5)=(6.-RSTAR)*(RSTAR**2.)*ESTAR3
ENDIF
ENDIF
c write(*,*) 'GGGGGGGGGGGGGGG'
50 IF(NOCC(6) .EQ. 0) GOTO 60
c write(*,*) 'HHHHHHHHHHHHHHH'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR3 .EQ. 0.))THEN
PSTOR(J,6)=0.
ELSE
DPSTOR=3.*LOG10(RSTAR)+LOG10(ESTAR3)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,6)=0.
ELSE
PSTOR(J,6)=(RSTAR**3.)*ESTAR3
ENDIF
ENDIF
c write(*,*) 'NNNNNNNNNNNNN'
60 IF(NOCC(7) .EQ. 0) GOTO 70
c write(*,*) 'IIIIIIIIIIIIIII'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR4 .EQ. 0.) .OR.
/ ((192.-144.*RSTAR+24.*(RSTAR**2.)-(RSTAR**3.))
/ .EQ. 0.))THEN
PSTOR(J,7)=0.
ELSE
DPSTOR=LOG10(ABS((192.-144.*RSTAR+24.*(RSTAR**2.)
/ -(RSTAR*3.))))+LOG10(RSTAR)+LOG10(ESTAR4)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,7)=0.
ELSE
PSTOR(J,7) =(192.-144.*RSTAR+24.*(RSTAR**2.)
/ -(RSTAR**3.))*RSTAR*ESTAR4
ENDIF
ENDIF
c write(*,*) 'OOOOOOOOOOOOOOOO'
70 IF(NOCC(8) .EQ. 0) GOTO 80
c write(*,*) 'JJJJJJJJJJJJJJJJJJ'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR4 .EQ. 0.) .OR.
/ ((80.-20.*RSTAR+(RSTAR**2.)) .EQ. 0.))THEN
PSTOR(J,8)=0.
ELSE
DPSTOR=LOG10(ABS(80.-20.*RSTAR+(RSTAR**2.)))
/ +2.*LOG10(RSTAR)+LOG10(ESTAR4)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,8)=0.
ELSE
PSTOR(J,8)=(80.-20.*RSTAR+(RSTAR**2.))*(RSTAR**2.)*ESTAR4
ENDIF
ENDIF
c write(*,*) 'PPPPPPPPPPPP'
80 IF(NOCC(9) .EQ. 0) GOTO 90
c write(*,*) 'KKKKKKKKKKKKKKKKKKK'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR4 .EQ. 0.) .OR.
/ ((12.-1.*RSTAR) .EQ. 0.)) THEN
PSTOR(J,9)=0.
ELSE
DPSTOR=LOG10(ABS(12.-1.*RSTAR))
/ +3.*LOG10(RSTAR)+LOG10(ESTAR4)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,9)=0.
ELSE
PSTOR(J,9)=(12.-1.*RSTAR)*(RSTAR**3.)*ESTAR4
ENDIF
ENDIF
90 IF(NOCC(10) .EQ. 0) GOTO 100
c write(*,*) 'LLLLLLLLLLLLLLLLLLLL'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR4 .EQ. 0.))THEN
PSTOR(J,10)=0.
ELSE
DPSTOR=4.*LOG10(RSTAR)+LOG10(ESTAR4)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,10)=0.
ELSE
PSTOR(J,10)=(RSTAR**4.)*ESTAR4
ENDIF
ENDIF
100 IF(NOCC(11) .EQ. 0) GOTO 110
C 100 IF(NOCC(11) .EQ. O) GOTO 115
c write(*,*) 'MMMMMMMMMMMMMMMMMM'
DPSTOR=0.
IF((RSTAR .EQ. 0.) .OR. (ESTAR5 .EQ. 0.) .OR.
/ ((9375.-7500.*RSTAR+1500.*(RSTAR*2.)
/ -100.*(RSTAR**3.)+2.*(RSTAR**4.)) .EQ. 0.))THEN
PSTOR(J,11)=0.
ELSE
DPSTOR=LOG10(ABS(9375.-7500.*RSTAR+1500.*(RSTAR**2.)
/ -100.*(RSTAR**3.)+2.*(RSTAR**4.)))
/ +LOG10(RSTAR)+LOG10(ESTAR5)
IF(DPSTOR .LT. ((-307.)+LOG10(4.3)))THEN
PSTOR(J,11)=0.
ELSE
PSTOR(J,11)=(9375.-7500.*RSTAR+1500.*(RSTAR**2.)
/ -100.*(RSTAR**3.)+2.*(RSTAR**4.))
/ *RSTAR*ESTAR5
ENDIF
ENDIF
C 115 WRITE(6,*)'J,PSTOR(J,1),PSTOR(J,2),PSTOR(J,3),PSTOR(J,4)',
C / J,PSTOR(J,1),PSTOR(J,2),PSTOR(J,3),PSTOR(J,4)
C WRITE(6.*)'J.PSTOR(J,5),PSTOR(J,6),PSTOR(J,7),PSTOR(J,8)',
C J,PSTOR(J,5),PSTOR(J,6),PSTOR(J,7),PSTOR(J,8)
C WRITE(6.*)'J.PSTOR(J,9),PSTOR(J,10),PSTOR(J.11)',
C / J,PSTOR(J.9),PSTOR(J,10),PSTOR(J,11)
110 CONTINUE
c stop
c write(*,*) 'QQQQQQQQQQQQQQQQQ'
DO 140 STATE=1,NSTATE
IF(NOCC(STATE) .EQ. 0) GOTO 140
NORM=0.
DO 120 J=1,NR
IF (ABS(PSTOR(J,STATE)) .LT. (0.656D-153)) GOTO 120
NORM=NORM+(PSTOR(J,STATE)**2)
120 CONTINUE
C write(*,*) 'NORM*DR',NORM*DR,norm
NORM=1./DSQRT(NORM*DR)
DO 130 J=1,NR
C DPLOOG=0.
DPLOOG=0.
IF((PSTOR(J,STATE) .EQ. 0.) .OR. (NORM .EQ. 0.))THEN
PSTOR(J,STATE)=0.
C write(*,*) 'if 1'
ELSE
C write(*,*) 'if 2-1'
C write(*,*) 'if 2-1',pstor(j,state),norm
DPLOOG=LOG10(ABS(PSTOR(J,STATE)))+LOG10(ABS(NORM))
C write(*,*) 'if 2-154654yjyt',pstor(j,state),norm
IF(DPLOOG .LT. ((-307.)+LOG10(4.3)))THEN
C write(*,*) 'if 3-1'
PSTOR(J,STATE)=0.
ELSE
C write(*,*) 'if 2-2'
PSTOR(J,STATE)=PSTOR(J,STATE)*NORM
ENDIF
C write(*,*) 'if 2'
ENDIF
C WRITE(6,*)'J,STATE,PSTOR(J,STATE)',J,STATE,PSTOR(J,STATE)
130 CONTINUE
140 CONTINUE
C WRITE(6,*)'RETURN Normalization'
RETURN
END
C*************************************************
SUBROUTINE ENERGY
DIMENSION PSTOR(0:1400,11),RHO(1400), NOCC(11),FOCK(1400,11),
/AML(11),PHI(0:1400),EXH(11)
COMMON /GENE1/DR,NR/GENE3/NSTATE/KIKAKU/NEE,PFLAG,NIT,VENTOT,
/VEETOT,VEXTOT,KTOT,ELAST,ETOT/KYOU1/ZE2,HBM/KYOU2/NOCC/KYOU3/AML
//JIGEN8/EXH/JIGEN2/RHO/JIGEN3/PHI/JIGEN4/FOCK/JIGEN7/ID
//JIGEN1/PSTOR
REAL*8 PSTOR,RHO,FOCK,DR,PHI,KTOT,VENTOT,VEETOT,VEXTOT,KEN,
/VEN,VEE,VEX,VCENT,PM,R,PZ,PZ2,HBM,ZE2,EXH,ELAST,ETOT
INTEGER*4 PFLAG,NR,NSTATE,NOCC,AML,NIT,LL1,NEE,STATE
CHARACTER*14 ID(11)
C WRITE(6,*) 'ENERGY'
C ***** calcurate density and Fock ****************
C WRITE(6,*)'call density'
CALL DENFO
C ******************************
C ****************************
C WRITE(6,*)'CALL Poisson '
CALL POISSON
C **********************************
C ****************************************
KTOT=0.
VENTOT=0.
VEETOT=0.
VEXTOT=0.
IF(PFLAG .EQ. 0)GOTO 20
WRITE(6,10)NIT
10 FORMAT(1H ,'---- Iteration -------- Round ',I5,' ----')
20 DO 70 STATE=1 ,NSTATE
C write(*,*) 'do loop70',state,Nstate,NOCC(state)
IF (NOCC(STATE) .EQ. 0) GOTO 70
KEN=0.
VEN=0.
VEE=0.
VEX=0.
VCENT=0.
LL1=AML(STATE)*(AML(STATE)+1)
C WRITE(*,*) LL1
PM=0.
DO 30 J=1,NR
R=DBLE(J)*DR
C write(*,*) 'do loop30',j,Nr,R
IF(ABS(PSTOR(J,STATE)) .LT. (0.656D-153)) THEN
PZ=0.
ELSE
PZ=PSTOR(J,STATE)
ENDIF
C write(*,*) 'PZ=',PZ,PM
C2022.1.23 PZ2=PZ**2.
PZ2=ABS(PZ)**2.
C2022.1.23 KEN=KEN+(PZ-PM)**2
KEN=KEN+ABS(PZ-PM)**2
VCENT=VCENT+PZ2*DBLE(LL1)/(R**2.)
VEN=VEN-PZ2/R
VEE=VEE+PHI(J)*PZ2
VEX=VEX+FOCK(J,STATE)*PZ
PM=PZ
30 CONTINUE
VCENT=VCENT*DR*HBM/2.
KEN=KEN*HBM/(2.*DR)
VEE=VEE*DR
VEN=VEN*ZE2*DR
VEX=VEX*DR
KEN=KEN+VCENT
EXH(STATE)=KEN+VEN+VEE+VEX
KTOT=KTOT+KEN*DBLE(NOCC(STATE))
VENTOT=VENTOT+VEN*DBLE(NOCC(STATE))
VEETOT=VEETOT+VEE*DBLE(NOCC(STATE))
VEXTOT=VEXTOT+VEX*DBLE(NOCC(STATE))
IF(PFLAG .EQ. 0)GOTO 70
WRITE(6,40)ID(STATE),NOCC(STATE),AML(STATE)
40 FORMAT(1H ,'orbit ',A14,' electron',I5,
/ ' arg. momentum',I5/)
WRITE(6,50)
50 FORMAT(4x,'momentum',4X,'a.force',6x,'r.force',6X,'e.force',
/ 11X,'V',11X,'e')
WRITE(6,60)KEN,VEN,VEE,VEX,VEN+VEE+VEX,EXH(STATE)
60 FORMAT(1H ,6E13.5/)
70 CONTINUE
VEETOT=VEETOT/2.
VEXTOT=VEXTOT/2.
ETOT=KTOT+VENTOT+VEETOT+VEXTOT
IF(PFLAG .EQ. 0) GOTO 120
WRITE(6,*)'*** Totales ***'
WRITE(6,80)
80 FORMAT(2X,'electron',2X,'momentum',5X,'a.force',7X,'r.force',7X,
/'e.force',12X,'V')
WRITE(6,90)NEE,KTOT,VENTOT,VEETOT,VEXTOT,VENTOT+VEETOT+VEXTOT
90 FORMAT(1H ,I6,2X,5E14.6/)
WRITE(6,100)ETOT
100 FORMAT(1H ,'Total energy = ',E15.8)
WRITE(6,110)ELAST
110 FORMAT(1H ,'Previous total energy = ',E15.8)
120 ELAST=ETOT
CC WRITE(6,*)'RETURN Energy'
RETURN
END
C*********Density and Fock****************
SUBROUTINE DENFO
DIMENSION RHO(1400),NOCC(11),PSTOR(0:1400,11),FOCK(1400,11),
/AML(11)
COMMON /GENE1/DR,NR/GENE2/E2/GENE3/NSTATE/MITUDO/MIX/COF/L1,L2,
/LAM,THREEJ/KYOU2/NOCC/KYOU3/AML/JIGEN1/PSTOR/JIGEN2/RHO
//JIGEN4/FOCK
REAL*8 RHO,MIX,PSTOR,FOCK,FAC,E2,SUM,R,RLAM,TERM,DR,DF,RLAM1,
/THREEJ,DFLOG,DTERM
INTEGER*4 NR,NSTATE,NOCC,AML,L1,L2,LAMSTART,LAMSTOP,LAM,STATE,
/STATE2
C WRITE(6,*)'density',log10(2.)
DO 10 J=1,NR
RHO(J)=(1.-MIX)*RHO(J)
10 CONTINUE
C WRITE(6,*)'densityytyutrjytmjj'
DO 90 STATE=1,NSTATE
IF (NOCC(STATE) .EQ. 0)GOTO 90
C WRITE(6,*)'densityyty not goto 90'
DO 20 J=1,NR
C write(*,*) 'PSTOR(J,STATE)',PSTOR(J,STATE),j,state
IF(ABS(PSTOR(J,STATE)) .LT. (0.656D-153))THEN
C WRITE(6,*)'rho(j)=rho(j)'
RHO(J)=RHO(J)
ELSE
C WRITE(6,*)'%%%%%%%%%%%%%%%%%%%%%%%%%%%'
RHO(J)=RHO(J)+MIX*DBLE(NOCC(STATE)*PSTOR(J,STATE)**2)
C WRITE(6,*)'%%%%%%%%%%%%%%%',j,nr
ENDIF
C WRITE(6,*)'J,STATE,RHO(J)',J,STATE,RHO(J)
FOCK(J,STATE)=0.
C WRITE(6,*)'J,STATE,RHO(J)'
20 CONTINUE
C WRITE(6,*)'STATE,DO 1 ',STATE
DO 80 STATE2=1,NSTATE
IF(NOCC(STATE2) .EQ. 0) GOTO 80
L1=AML(STATE)
L2=AML(STATE2)
LAMSTART=ABS(L1-L2)
LAMSTOP=L1+L2
DO 70 LAM=LAMSTART,LAMSTOP,2
C**************************
C *****************************
C WRITE(6,*)'CALL EXCHANGE LAM=',LAM
CALL KOUKAN
FAC=-E2/2.*DBLE(NOCC(STATE2))*THREEJ
SUM=0.
DO 40 I=1,NR
DTERM=0.
DFLOG=0.
R=DBLE(I)*DR
RLAM=R**DBLE(LAM)
C WRITE(*,*) 'LAM=',LAM,'i=',i
IF((PSTOR(I,STATE2) .EQ. 0.) .OR.
/ (PSTOR(I,STATE) .EQ. 0.)) THEN
TERM=0.
C write(*,*)'density2-0 if1'
ELSE
C write(*,*)'density2-0 else1',I,STATE
DTERM=LOG10(ABS(PSTOR(I,STATE2)))
/+LOG10(ABS(PSTOR(I,STATE)))
/+LOG10(RLAM)-LOG10(2.)
IF(DTERM .LT. ((-307.)+LOG10(4.3)))THEN
C write(*,*)'density2-0 if2'
TERM=0.
ELSE
c write(*,*)'density2-0 else2'
TERM=PSTOR(I,STATE2)*PSTOR(I,STATE)*RLAM/2.
ENDIF
ENDIF
c write(*,*)'density2'
SUM=SUM+TERM
c write(*,*)'density2-1'
IF((PSTOR(I,STATE2) .EQ. 0.) .OR. (SUM .EQ. 0.)
/ .OR. (FAC .EQ. 0.)) GOTO 30
DFLOG=LOG10(ABS(PSTOR(I,STATE2)))+LOG10(ABS(SUM))
/ +LOG10(ABS(FAC))+LOG10(DR)
/ -LOG10(RLAM)-LOG10(R)
c write(*,*)'density3'
IF(DFLOG .LT. ((-307.)+LOG10(4.3))) GOTO 30
DF=PSTOR(I,STATE2)*FAC*SUM*DR/(RLAM*R)
FOCK(I,STATE)=FOCK(I,STATE)+DF
30 SUM=SUM+TERM
40 CONTINUE
SUM=0.
C write(*,*)'density4'
DO 60 I=1,NR
II=NR+1-I
DTERM=0.
DFLOG=0.
R=DBLE(11)*DR
RLAM1=R**DBLE((LAM+1))
IF((PSTOR(II,STATE2) .EQ. 0.) .OR.
/ (PSTOR(II,STATE) .EQ. 0.)) THEN
TERM=0.
ELSE
DTERM=LOG10(ABS(PSTOR(II,STATE2)))
/ +LOG10(ABS(PSTOR(II,STATE)))
/ -LOG10(RLAM1)-LOG10(2.)
IF(DTERM .LT. ((-307.)+LOG10(4.3))) THEN
TERM=0.
ELSE
TERM=PSTOR(II,STATE2 )*PSTOR(II,STATE)/RLAM1/2.
ENDIF
ENDIF
SUM=SUM+TERM
IF((PSTOR(II,STATE2) .EQ. 0.) .OR. (SUM .EQ. 0.)
/ .OR. (FAC .EQ. 0.)) GOTO 50
DFLOG=LOG10(ABS(PSTOR(II,STATE2)))+LOG10(ABS(SUM))
/ +LOG10(ABS(FAC))+LOG10(DR)
/ +LOG10(RLAM1)-LOG10(R)
IF(DFLOG .LT. ((-307.)+LOG10(4.3))) GOTO 50
DF=PSTOR(II,STATE2)*FAC*SUM*DR*RLAM1/R
FOCK(II,STATE)=FOCK(II,STATE)+DF
50 SUM=SUM+TERM
60 CONTINUE
70 CONTINUE
80 CONTINUE
90 CONTINUE
c WRITE(6,*)'RETURN Density'
RETURN
END
- 前ページ
- 次ページ