Advertisement
senj0h

hfock source needs debugging

Jan 7th, 2016
152
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. C     THE ATOM HARTREE FOCK
  2. C 1-3-2016 Eric Lovejoy removes impicit declaration
  3. C     organized code, and fixed typos
  4.       PROGRAM HFOCK
  5.       implicit none
  6.  
  7.       INTEGER I
  8.       REAL A, B, C, E, F, T
  9.       REAL C1, C2, ET
  10.       REAL F11,F12,F22,H11,H12,H22,R11,R12,R22
  11.       REAL S12,T11,T12,T22,V11,V12,V22
  12.       REAL DIFF,EOLD,PLUS,PROD,C1C2
  13.       REAL ZETA1,ZETA2
  14.       REAL ER1111,ER1112,ER1122,ER1212,ER1222,ER2222
  15.  
  16.       ZETA1 = 2.173171
  17.       ZETA2 = 1.188530
  18.       WRITE (*, '( "Zeta1:", F7.4 / "Zeta2:", F7.4 )' )  ZETA1, ZETA2
  19.  
  20. C COMPUTE INTEGRALS
  21.       PLUS=ZETA1+ZETA2
  22.       PROD=ZETA1*ZETA2
  23.       DIFF=ZETA1-ZETA2
  24.       S12=8.*PROD**1.5/PLUS**3
  25.       T11=0.5*ZETA1**2
  26.       T22=0.5*ZETA2**2
  27.       T12=4.*PROD**2.5/PLUS**3
  28.       V11=-2.*ZETA1
  29.       V22=-2.*ZETA2
  30.       V12=-8.*PROD**1.5/PLUS**2
  31.       H11=T11+V11
  32.       H12=T12+V12
  33.       H22=T22+V22
  34.       ER1111=5.*ZETA1/8.
  35.       ER2222=5.*ZETA2/8.
  36.       T=(DIFF/PLUS)**2
  37.       ER1122=PLUS*(1.-T)*(5.-T)/16.
  38.       ER1212=4.*PROD**3/PLUS**5
  39.       F=0.25*PROD**1.5/PLUS**3
  40.       T=(DIFF/(3.*ZETA1+ZETA2)**2)
  41.       ER1112=F*(3.*ZETA1+ZETA2)*(1.-T)*(5.-T)
  42.       T=(DIFF/(3.*ZETA2+ZETA1)**2)
  43.       ER1222=F*(3.*ZETA2+ZETA1)*(1.-T)*(5.-T)
  44.       R11=1.
  45.       R12=0.
  46.       R22=0.
  47.       EOLD=0
  48. c variables
  49. c S12 = 0.874115, T11 = 2.361336, T12 = 1.128867, T22 = 0.706302
  50. c V11 = -4.346342,      V12 = -2.938515,    V22 = -2.377060
  51. c ER1111, = 1.358232,   ER1222  = 0.756250
  52. c ER1112  = 1.032210        ER2222  = 0.742831
  53. c ER1122, = 0.943928        ER1212  = 0.802688
  54.       WRITE(*,*) '***********************************************'
  55.       WRITE(*,*) '*****System Variables: Actual vs Expected*****'
  56.       WRITE(*,*) '***********************************************'  
  57.       WRITE(*,*) 'S12',S12, ' Expected:0.874115'
  58.       WRITE(*,*) 'T11',T11, ' Expected:2.361336'
  59.       WRITE(*,*) 'T12',T12, ' Expected:1.128867'
  60.       WRITE(*,*) 'T22',T22, ' Expected:0.706302'
  61.       WRITE(*,*) 'V11',V11, ' Expected:-4.346342'
  62.       WRITE(*,*) 'V12',V12, ' Expected:-2.938515'
  63.       WRITE(*,*) 'V22',V22, ' Expected:-2.377060'
  64.       WRITE(*,*) 'ER1111',ER1111, ' Expected:1.358232'
  65.       WRITE(*,*) 'ER1222',ER1222, ' Expected:0.756250'
  66.       WRITE(*,*) 'ER1112',ER1112, ' Expected:1.032210'
  67.       WRITE(*,*) 'ER2222',ER2222, ' Expected:0.742831'
  68.       WRITE(*,*) 'ER1122',ER1122, ' Expected:0.943928'
  69.       WRITE(*,*) 'ER1212',ER1212, ' Expected:0.802688'
  70.  
  71.       WRITE(*,*) ' '
  72.  
  73.       WRITE(*,*) '***********************************************'
  74.       WRITE(*,*) '* ITERATION ORBITAL ENERGY TOTAL E EIGENVECTOR*'
  75.       WRITE(*,*) '***********************************************'
  76. C  START HARTREE-FOCK LOOP
  77.       DO 3 I=1,15
  78.       F11=H11+R11*ER1111+2.*R12*ER1112+R22*(2.*ER1122-ER1212)
  79.       F12=H12+R11*ER1112+R12*(3.*ER1212-ER1122)+R22*ER1222
  80.       F22=H22+R11*(2.*ER1122-ER1212)+2.*R12*ER1222+R22*ER2222
  81.       A=1.-S12**2
  82.       B=2.*F12*S12-F11-F22
  83.       C=F11*F22-F12**2
  84.       E=(-B-(B**2-4.*A*C)**0.5)/(2.*A)
  85.       C1C2=-(F12-E*S12)/(F11-E)
  86.       C2=1./(C1C2**2+2.*C1C2*S12+1.)**0.5
  87.       C2=1./(C1C2**2+2.*C1C2*S12+1.)**0.5
  88.       C1=C1C2*C2
  89.       R11=C1**2
  90.       R12=C1*C2
  91.       R22=C2**2
  92.       ET=E+R11*H11+2.*R12+R22*H22
  93.       WRITE(*,*) I,E,ET,C1,C2
  94.       IF((ABS(ET-EOLD)).LT.1.E-6) STOP
  95.       EOLD=ET
  96.  3    CONTINUE
  97.       STOP
  98.       END PROGRAM HFOCK
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement