API tools faq
paste
Advertisement
j0h

integration.f

j0h
Jun 20th, 2024
537
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 61.21 KB | None | 0 0
raw download clone embed print report
  1. C   Copy of INTEGRAL FORTRAN WRITTEN BY Jay Rasaiah
  2. C   MODIFIED INPUT TO RUN ON IMAC by E.LOVEJOY June 18 2024, Rev 8.01  
  3. C
  4. C   Compiled with:
  5. C   gfortran -g -fcheck=all -std=legacy  -o integral integral.f
  6. C   -fcheck=all enables array bound checks at run time
  7. C **********************************************************************INT00020
  8. C  THIS IS VERSION 1.57 WHICH USES LESS MEMORY AND IS FASTER THAN      *INT00030
  9. C  VERSION 1.55 AND CALC THE SHORT-RANGED DIR CORR FUNCTION(DCS)        INT00040
  10. C  PROGRAM RESTRICTED TO SYMMETRICALLY CHARGED ELECTROLYTES.           *INT00050
  11. C  SOLVES HNC EQ BY ITERATIONUSING ISEINGS FAST FOURIER TRANSFORM FFT3S*INT00060
  12. C  IF (CONC.EQ.0) READS A NEW SET OF ION AND WELL SIZE PARAMETERS.     *INT00070
  13. C  IF REM=0.0, NO REMAINDER ADDED TO G(I,J), T(2,4), G4(I,J)           *INT00080
  14. C  IF REM.NE.0.0, DEBYE-HUCKEL REMAINDER ADDED                         *INT00090
  15. C  PRIMITIVE MODEL PLUS TWO SQUARE WELLS                               *INT00100
  16. C  SPAC IS NEAREST FRACTION OF A(1,2) THAT IS CLOSEST TO               *INT00110
  17. C  MAX(COEFF1*DEBYE/FNN,COEFF2*A(1,2)/FNN)                             *INT00120
  18. C  DELA = DLN(ACTCO)/DLNC, CALC FROM COMPRESSIBILITY EQUATION.         *INT00130
  19. C  DELB = LIMIT(DELB) AS BETA GOES TO ZERO                             *INT00140
  20. C  BETA MEAS DEVIATION FROM ELECTRONEUTRALITY                          *INT00150
  21. C  OSM IS CALCULATED FROM THE PRESSURE EQUATION                        *INT00160
  22. C  XENERGY AND DFX/DP FROM INDEP EQNS ASSUMING NO TEMP AND PRESS DEP   *INT00170
  23. C  OF THE SHORT RANGE POTENTIAL XENERGY IS IN CALS PER MOLE OF SOLUTE  *INT00180
  24. C  XVOLUME IN CC PER MOLE OF SOLUTE                                    *INT00190
  25. C  DRF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE HNC EQN             *INT00200
  26. C  DRDF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE DHLL               *INT00210
  27. C  CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12-C(1)+C(2)                   *INT00220
  28. C  CONVGB = 100*(G(1,2)-G(2,1))/GM12                                   *INT00230
  29. C  DCS(=ZETA) IS THE SHORT-RANGED DIRECT CORRELATION FUNCTION          *INT00240
  30. C  CONVERGENCE CRITERION USEFUL ONLY FOR UNSYMMETRICAL SYSTEMS         *INT00250
  31. C  IF EX.EQ.2.0, PRINTS XR,DCS,H                                       *INT00260
  32. C  IF EX.EQ.3.0, PRINTS XR,DCS,H AND CHARGE DENSITY                    *INT00270
  33. C  IF EX.EQ.0.0, QUITS                                                 *INT00280
  34. C  IF CHARGE = 0.0, REPEATS CALC WITH SAME SPACING AS CHARGED SYSTEM   *INT00290
  35. C  BUT WITH Z1 = Z2 = 0.0.FICTITIOUS CHARGE DENSITY ALSO CALC          *INT00300
  36. C  IF ICC.GT.0, ITERATION BEGINS WITH TAU OF PRECEDING CONC            *INT00310
  37. C  IF ICC.LT.0, READ TAU GENERATED IN A PREVIOUS EXECUTION             *INT00320
  38. C  IF ICC.EQ.0, TAU CAL IN LAMBDA APPROXIMATION                        *INT00330
  39. C  WRITES TAU ON NDSK2 IN ALL CASES                                    *INT00340
  40. C  SIX SUBRUTNIES ADJ,TRAP,SIMPSON,FFT3S,SETFT,RSA
  41. C  NIONS is number positive and negative ions 2 for NaCl and ZnSO4
  42.  
  43. C **************************DIMENSIONS**********************************INT00350
  44.       INTEGER MNM
  45.       PARAMETER (MNM = 2048)
  46.       REAL*8 DATE,TIME                                                  INT00360
  47.       DIMENSION A(2,2)
  48.       DIMENSION A1(2,2)
  49.       DIMENSION B(80)
  50.       DIMENSION BA(2,2)
  51.       DIMENSION B1(2,2)
  52.       DIMENSION C(2)
  53.       DIMENSION CB(2,2)
  54.       DIMENSION CONV(MNM,2,2)
  55.       DIMENSION C1(2,2)
  56.       DIMENSION DC(MNM,2,2)
  57.       DIMENSION DRF(MNM)
  58.       DIMENSION DRDF(MNM)
  59.       DIMENSION E(2,2)
  60.       DIMENSION FNU(2)
  61.       DIMENSION F2(2)
  62.       DIMENSION F4(2)
  63.       DIMENSION G(2,2)
  64.       DIMENSION GA(2,2)
  65.       DIMENSION GB(2,2)
  66.       DIMENSION GC(2,2)
  67.       DIMENSION GM(2,2)
  68.       DIMENSION GW(2,2)
  69.       DIMENSION G4(2,2)
  70. C      DIMENSION H(MNM,2,2)
  71.       REAL H(MNM,2,2)
  72.       DIMENSION LINT(2,2)
  73.       DIMENSION MINT(2,2)
  74.       DIMENSION NINT(2,2)      
  75.       DIMENSION PHI(MNM,2,2)
  76.       DIMENSION P1(MNM,2,2)  
  77.       DIMENSION Q(MNM,2,2)
  78.       DIMENSION Q1(MNM,2,2)
  79.       DIMENSION RHO(2,MNM)
  80.       DIMENSION RQO(2,MNM)
  81.       DIMENSION R2(2)
  82.       DIMENSION R4(2)
  83.       DIMENSION T(2,2)
  84.       DIMENSION TAU(MNM,2,2)
  85.       DIMENSION TCONV(MNM,2,2)
  86.       DIMENSION TH(MNM,2,2)
  87.       DIMENSION TM(2,2)
  88.       DIMENSION TPHI(MNM,2,2)
  89.       DIMENSION TQ(MNM,2,2)
  90.       DIMENSION TTAU(MNM,2,2)
  91.       DIMENSION TW(2,2)
  92.       DIMENSION TZETA(MNM,2,2)
  93.       DIMENSION W(2,2)
  94. C      DIMENSION XR(MNM)
  95.       REAL XR(MNM)
  96.       DIMENSION Y(2)
  97.       DIMENSION Y1(MNM)
  98.       DIMENSION Z(2)
  99.       DIMENSION ZA(2)
  100.       DIMENSION ZETA(MNM,2,2)
  101.    
  102.       EQUIVALENCE (H,TH,PHI,TPHI,TTAU)                                  
  103.       EQUIVALENCE (ZETA,TZETA,RHO,DRF,DC)                              
  104.       EQUIVALENCE (TAU,P1,CONV,TCONV)                                  
  105.       EQUIVALENCE (RQO,DRDF)                                            
  106.       EQUIVALENCE (XR,Y1)                                              
  107.       EQUIVALENCE (B(6),Z(1)),(B(8),D),(B(9),TEMP),(B(10),A1(1,1)),    
  108.      *(B(14),COEFF1),(B(15),DEBYE),(B(16),SPAC),(B(17),FNN),(B(18),CONC)
  109.      *,(B(19),DELA),(B(20),OSM),(B(21),CONTAC),(B(22),REMAIN),          
  110.      *(B(23),FITER),(B(24),G(1,1)),(B(28),T(1,1)),(B(32),FNGARB),      
  111.      *(B(33),DELB),(B(34),DIFF),(B(35),CONVG1),(B(36),CONVG2),          
  112.      *(B(37),GA(1,1)),(B(5),DT),(B(4),DP),(B(41),GB(1,1)),              
  113.      *(B(45),G4(1,1)),(B(49),F2(1)),(B(51),F4(1)),(B(52),COEFF2)        
  114.       COMMON /JAY/ H,XR
  115. C******************************* Hunting for issues********************
  116. C     Gfortran has some other nifty flags for ID of various issues:
  117. C     -ffpe-trap=invalid,zero,overflow
  118. C     gfortran -g -ffpe-trap=invalid,zero,overflow -std=legacy -fcheck=all
  119.  
  120. C******************************* DATE TIME FIX ************************
  121.       CHARACTER*25 date_string
  122.       CHARACTER*100 command
  123.       command = 'date "+%Y-%m-%d %H:%M:%S" > date.txt'
  124.       CALL SYSTEM(command)
  125.       OPEN(UNIT=1, FILE='date.txt', STATUS='OLD')
  126.       READ(1, '(A)') date_string
  127.       CLOSE(1)
  128.  
  129.  
  130.  
  131. C ******************************* CONSTANTS ****************************INT00600
  132. C     PROGRAM CONSTANTS                                                *INT00610
  133.       DATA PI/3.141593/                                                 INT00620
  134.       PRINT *,  date_string
  135. C     CALL TIMER(DATE,TIME,IVIRT,ITOT)                                  INT00630
  136.       NRDR  = 5                                                         INT00640
  137.       NPRT  = 7                                                         INT00650
  138.       NPUN  = 7                                                         INT00660
  139.       NDSK  = 8                                                         INT00670
  140.       NDSK2 = 9                                                         INT00680
  141.       NIONS = 2                                                         INT00690
  142.       PI2   = 2.*PI                                                     INT00700
  143.       PI4   = 4.*PI                                                     INT00710
  144.       PI43  = PI4/3.                                                    INT00720
  145.       PIPI  = PI*PI                                                     INT00730
  146. C **********************************************************************INT00740
  147. C     CONSTANTS FOR SUBROUTINE FFT3S                                   *INT01430
  148.       F0 = 0.                                                           INT01440
  149.       RL = SPAC1*FNN                                                    INT01450
  150. C **********************************************************************INT01460
  151.  
  152.      
  153. C*****used to stop crashing of  H MINT PHI TPHI TCONV CONV ZETA*******
  154.       I = 1
  155.       J = 1
  156.       IF (J .LT. 1) THEN
  157.       J = J + 1
  158.       END IF
  159.       IF (I .LT. 1) THEN
  160.       I = I + 1
  161.       END IF
  162.      
  163.       DO 10 I = 1,80                                                    INT00750
  164. 10    B(I) = 0.0                                                        INT00760
  165.       B(3) = 1.57                                                       INT00770
  166. C20   READ(NRDR,9001) M1                                                INT00780
  167. C20   READ(NRDR,*) M1                                                   INT00780
  168. C     READ(NRDR,*) (Z(I),FNU(I),I=1,NIONS)
  169.       M1 = 11        
  170. C      Z(1) = 1.0
  171.       Z(1) = 2.0
  172.       Z(2) = -1.0                                                       INT00790
  173.       ZE = Z(1)                                                         INT00800
  174.       ZA(1) = Z(1)                                                      INT00810
  175.       ZA(2) = Z(2)                                                      INT00820
  176.       FNU(1) = 1.0
  177.       FNU(2) = 2.0
  178. C25   READ(NRDR,*) D,TEMP,DT,DP                                        
  179.        D = 78.358
  180.        T = 298.15
  181.        TEMP = 298.15
  182.        DT = 0.0
  183.        DP = 0.0
  184. C30   READ(NRDR,*) A(1,1),BA(1,1),CB(1,1),A(1,2),BA(1,2),CB(1,2),A(2,   INT00840
  185. C    *2),BA(2,2),CB(2,2)                                                INT00850      
  186.       A(1,1) = 4.2
  187.       A(1,2) = 4.2
  188.       A(2,1) = A(1,2)                                                   INT00860
  189.       A(2,2) = A (1,1)
  190.       BA(1,1) = 0
  191.       BA(2,2) = BA(1,1)
  192.       BA(1,2) = 0.0
  193.       BA(2,1) = BA(1,2)
  194.       CB(1,1) = 0.0
  195.       CB(2,2) = CB(1,1)
  196.       CB(1,2) = 0.0                                                     INT00870
  197.       CB(2,1) = CB(1,2)                                                 INT00880
  198. C35   READ(NRDR,*) E(1,1),W(1,1),E(1,2),W(1,2),E(2,2),W(2,2)            INT00890
  199.       E(1,1)= 0.0
  200.       E(1,2)= 0.0
  201.       E(2,1) = E(1,2)                                                   INT00900
  202.       E(2,2) = E(1,1)                                                   INT00900
  203.       W(1,1) = 0.0
  204.       W(1,2) = 0.0
  205.       W(2,1) = W(1,2)
  206.       W(2,2) = W(1,1)
  207.  
  208. C40   READ(NRDR,*)CONC,COEFF1,COEFF2,EX,REM,CHARGE,MAX,I2,NCD,IC,       INT00920
  209. C    *NR,IR,ICC
  210.       CONC = 0.1                                                        INT00930
  211.       COEFF1=1.0
  212.       COEFF2=1.0
  213.       EX = 3.0
  214.       REM = 1.0
  215.       CHARGE = 1.0
  216.       MAX = 6
  217.       I2 =2
  218.       ICC= 1
  219. C     NCD =
  220. C     IC =
  221.       IR = 1
  222.       NR = 5
  223. C     REWIND NDSK2                                                      INT00940
  224. C     IF(ICC.LT.0)READ(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=INT00950
  225. C    *1,N)                                                              INT00960
  226. C     REWIND NDSK2                                                      INT00970
  227.       IF (EX.EQ.0.0) GO TO 5050                                         INT00980
  228. C     IF(CONC.EQ.0.0) GO TO 30                                          INT00990
  229. 9001  FORMAT(I3)                                                        INT01000
  230. 9002  FORMAT(8F10.5)                                                    INT01010
  231. 9003  FORMAT(4F10.5)                                                    INT01020
  232. 9004  FORMAT(9F5.2)                                                     INT01030
  233. 9005  FORMAT(6F10.5)                                                    INT01040
  234. 9006  FORMAT(1F10.7,5F5.2,7I5)                                          INT01050
  235.       WRITE(NPRT,9100)B(3), date_string
  236. 9100  FORMAT(1H1, 38X, 'HNC VERSION ', F4.2, ' - IBM 370/148 -  ', A25)
  237. C      WRITE(NPRT,9100)B(3),DATE,TIME                                    INT01060
  238. C9100  FORMAT(1H1,38X,'HNC VERSION ',F4.2,' - IBM 370/148 -  ',A8,2X,A8) INT01070
  239. C**** N=2**M1=2048 for M1=11, ***************************************
  240.       N=2**M1                                                           INT01090
  241.       N2=N/2                                                            INT01100
  242.       FNN  = N                                                          INT01120
  243.       FC   = 0.0                                                        INT01130
  244.       FNO  = 0.0                                                        INT01140
  245.       CC   = 0.0                                                        INT01150
  246.       DO 50 I = 1,NIONS                                                 INT01160
  247.       FNO  = FNU(I)+FNO                                                 INT01170
  248.       C(I) = FNU(I)*CONC*6.0233*1.E-04                                  INT01180
  249.       CC   = C(I)+CC                                                    INT01190
  250. 50    FC   = C(I)*Z(I)**2+FC
  251. C********************************************************************** INT01200
  252. C     Electron charge e=4.8032x10E10 cm**3/2 g**1/2 s**-2 (esu)
  253. C     Boltzconst k=1.3803x10E-16 cm**2 s**2/K                           INT01200
  254. C     Bjerrum length |z+z+|e**2/DkT=16.714x1.E+04/DT
  255. C***********************************************************************
  256.       FL1  = 16.7094*1.E+04/(D*TEMP)                                    INT01210
  257.       FL   = PI4*FL1                                                    INT01220
  258.       FK   = (FC*FL)**0.5                                               INT01230
  259.       FKFK  = FK*FK                                                     INT01240
  260.       DEBYE = 1./FK                                                     INT01250
  261.       FKA  = A(1,2)/DEBYE                                               INT01260
  262. C     REWIND NDSK                                                       INT01270
  263. C     WRITE(NDSK,9101)A(1,1),A(1,2),A(2,2),D,TEMP,CONC,FKA,DEBYE,Z(1),Z(INT01280
  264. C    *2)                                                                INT01290
  265. 9101  FORMAT(3F8.5,2F8.4,3F10.6,2F5.1)                                  INT01300
  266.       SPAC1 = COEFF1*DEBYE/FNN                                          INT01310
  267.       SPAC2 = COEFF2*A(1,2)/FNN                                         INT01320
  268.       IF(SPAC2.GT.SPAC1) SPAC1=SPAC2                                    INT01330
  269.       S1 = A(1,2)/SPAC1                                                 INT01340
  270.       L = S1+1.                                                         INT01350
  271.       S2 = L                                                            INT01360
  272.       IF (S2-S1.GT.0.5)  GO TO 60                                       INT01370
  273.       SPAC1 = A(1,2)/S2                                                 INT01380
  274.       GO TO 70                                                          INT01390
  275. 60    SPAC1 = A(1,2)/(S2-1.)                                            INT01400
  276. 70    CONTINUE                                                          INT01410
  277. C **********************************************************************INT01420
  278. C     CONSTANTS FOR SUBROUTINE FFT3S                                   *INT01430
  279.       F0 = 0.                                                           INT01440
  280.       RL = SPAC1*FNN                                                    INT01450
  281. C     FX,FY and FZ are Debye Huckel remainders for ?
  282. C **********************************************************************INT01460
  283.       FI = 1.0                                                          INT01470
  284.       SPAC = SPAC1                                                      INT01480
  285.       DO 80 K = 1,N                                                     INT01490
  286.       XR(K) = SPAC1*FI                                                  INT01500
  287. 80    FI = FI+1.                                                        INT01510
  288.       XN = XR(N)                                                        INT01520
  289.       FKN = FK*XN                                                       INT01530
  290.       IF (REM.EQ.0.0) GO TO 90                                          INT01540
  291.       FX = FL*(1.+FKN)*EXP(-FKN)/FKFK                                   INT01550
  292.       FY = FL*EXP(-FKN)/FK                                              INT01560
  293.       FZ = FL*EXP(-FKN)*(FKN**3+3.*FKN**2+6.*FKN+6.)/(FKFK*FKFK)        INT01570
  294.       GO TO 100                                                         INT01580
  295. 90    FX = 0.0                                                          INT01590
  296.       FY = 0.0                                                          INT01600
  297.       FZ = 0.0                                                          INT01610
  298. 100   CONTINUE                                                          INT01620
  299.       DO 110 I = 1,NIONS                                                INT01630
  300.       DO 110 J = 1,NIONS                                                INT01640
  301.       CALL ADJ (A(I,J),A1(I,J),MINT(I,J),SPAC1)                         INT01650
  302.       CALL ADJ (BA(I,J),B1(I,J),NINT(I,J),SPAC1)                        INT01660
  303.       CALL ADJ (CB(I,J),C1(I,J),LINT(I,J),SPAC1)                        INT01670
  304. 110   CONTINUE                                                          INT01680
  305.       GO TO 130                                                         INT01690
  306. 120   PRINT *, date_string
  307.  
  308. C 120  CALL TIMER(DATE,TIME,IVIRT,ITOT)                                 INT01700
  309.       WRITE(NPRT,9100)B(3),date_string                                  INT01710
  310. C      WRITE(NPRT,9100)B(3),DATE,TIME                                    INT01710
  311.       Z(1)   =  0.0                                                     INT01720
  312.       Z(2)   =  0.0                                                     INT01730
  313.       ZE     =  0.0                                                     INT01740
  314.       CHARGE = 20.0                                                     INT01750
  315.       REM    =  0.0                                                     INT01760
  316.       DEBYE  =  0.0                                                     INT01770
  317.       BETA   =  0.0                                                     INT01780
  318.       CONVGA =  0.0                                                     INT01790
  319. 130   WRITE(NPRT,9103) Z(1),Z(2),D,TEMP,DT,DP,A1(1,1),A1(2,2),A1(1,2),  INT01800
  320.      *B1(1,1),B1(2,2),B1(1,2),C1(1,1),C1(2,2),C1(1,2),E(1,1),E(2,2),    INT01810
  321.      *E(1,2),W(1,1),W(2,2),W(1,2),COEFF1,COEFF2,SPAC1,N,REM,CONC,DEBYE  INT01820
  322. 9103  FORMAT('Z1=',F5.1,3X,'Z2=',F5.1,/,'DIEL CONST=',F7.3,3X,'TEMP=',  INT01830
  323.      *F7.3,'DLND/DLNT=',F10.5,3X,'DLND/DP=',F7.3,/,'A11=',F6.3,3X,'A22='INT01840
  324.      *,F6.3,3X,'A12=',F6.3,/,'B11=',F6.3,3X,'B22=',F6.3,3X,'B12=',F6.3,/INT01850
  325.      *'C11=',F6.3,3X,'C22=',F6.3,3X,'C12=',F6.3,/,                      INT01860
  326.      *'E11=',F6.3,3X,'E22=',F6.3,3X,'E12=',F6.3,/,                      INT01870
  327.      *'W11=',F6.3,3X,'W22=',F6.3,3X,'W12=',F6.3,/,                      INT01880
  328.      *'COEFF1=',F7.4,3X,'COEFF2=',F7.4,3X,'SPAC1=',F7.4,3X,'N=',I4,3X,  INT01890
  329.      *'REM=',F7.4,/,'CONC=',F10.6,3X,'DEBYE=',F10.5,/)                  INT01900
  330. C *** CALC OF Q,Q1,AND TQ                                               INT01910
  331. C     WRITE(6,5558)                                                     INT01920
  332. 5558  FORMAT(1X,'HO')                                                   INT01930
  333.       DO 180 I = 1,NIONS                                                INT01940
  334.       DO 180 J = 1,NIONS                                                INT01950
  335.       AB = A1(I,J)+B1(I,J)                                              INT01960
  336.       BC = AB+C1(I,J)                                                   INT01970
  337.       DO 180 K = 1,N                                                    INT01980
  338.       Q(K,I,J) = -Z(I)*Z(J)*FL*EXP(-FK*XR(K))/(PI4*XR(K))               INT01990
  339.       IF (XR(K) .LT. A1(I,J)) GO TO 170                                 INT02000
  340.       IF(ABS(XR(K)-AB).LT.1.E-04) GO TO 160                             INT02010
  341.       IF ((XR(K)-AB).LT.1.E-07) GO TO 150                               INT02020
  342.       IF(ABS(XR(K)-BC).LT.1.E-04) GO TO 140                             INT02030
  343.       IF((XR(K)-BC).LT.1.E-07) GO TO 160                                INT02040
  344. 140   Q1(K,I,J) = EXP(Q(K,I,J))                                         INT02050
  345.       GO TO 180                                                         INT02060
  346. 150   Q1(K,I,J) = EXP(Q(K,I,J)-E(I,J))                                  INT02070
  347.       GO TO 180                                                         INT02080
  348. 160   Q1(K,I,J) = EXP(Q(K,I,J)-W(I,J))                                  INT02090
  349.       GO TO 180                                                         INT02100
  350. 170   Q1(K,I,J) = 0.0                                                   INT02110
  351. 180   CONTINUE                                                          INT02120
  352. C     WRITE(6,5557)                                                     INT02130
  353. 5557  FORMAT(1X,'HA')                                                   INT02140
  354.       DO 190 I = 1,NIONS                                                INT02150
  355.       DO 190 J = 1,NIONS                                                INT02160
  356. 190   CALL FFT3S(F0,Q(1,I,J),M1,RL,TQ(1,I,J),1)                         INT02170
  357. C***************** Distributution ************************************
  358. C     CALC OF APPROX DISTRIBUTION FUNCTIONS after counter I1 =0        *INT02190
  359.  
  360. 1010  I1 = 0                                                            INT02210
  361.       IF (ICC.NE.0) GO TO 1070                                          INT02220
  362.       DO 1020 I = 1,NIONS                                               INT02230
  363.       DO 1020 J = 1,NIONS                                               INT02240
  364.       DO 1020 K = 1,N                                                   INT02250
  365. 1020  PHI(K,I,J) = Q1(K,I,J)-1.                                         INT02260
  366.       DO 1030 I = 1,NIONS                                               INT02270
  367.       DO 1030 J = 1,NIONS                                               INT02280
  368.       PHI(MINT(I,J),I,J) = (PHI(MINT(I,J),I,J)-1.)*.5                   INT02290
  369. 1030  CALL FFT3S(F0,PHI(1,I,J),M1,RL,TPHI(1,I,J),1)                     INT02300
  370.       DO 1040 I = 1,NIONS                                               INT02310
  371.       DO 1040 J = 1,NIONS                                               INT02320
  372.       DO 1040 K = 1,N                                                   INT02330
  373.       TCONV(K,I,J) = 0.0                                                INT02340
  374.       DO 1040 L = 1,NIONS                                               INT02350
  375. 1040  TCONV(K,I,J) = TPHI(K,I,L)*TPHI(K,L,J)*C(L)+TCONV(K,I,J)          INT02360
  376.       DO 1050 I = 1,NIONS                                               INT02370
  377.       DO 1050 J = 1,NIONS                                               INT02380
  378. 1050  CALL FFT3S(F0,TCONV(1,I,J),M1,RL,CONV(1,I,J),-1)                  INT02390
  379.       DO 1060 I = 1,NIONS                                               INT02400
  380.       DO 1060 J = 1,NIONS                                               INT02410
  381.       DO 1060 K = 1,N                                                   INT02420
  382. 1060  TAU(K,I,J) = CONV(K,I,J)+Q(K,I,J)*XR(K)*FK*.5                     INT02430
  383. 1070  DO 1080 I = 1,NIONS                                               INT02440
  384.       DO 1080 J = 1,NIONS                                               INT02450
  385.       DO 1075 K = 1,N                                                   INT02460
  386. 1075  H(K,I,J) = Q1(K,I,J)*(1.+TAU(K,I,J))-1.                           INT02470
  387.       IF (MINT(I,J) < 1) THEN
  388.       MINT(I,J) = 1
  389.       END IF              
  390.       H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5                       INT02480
  391.  
  392.       DO 1080 K = 1,N                                                   INT02490
  393. 1080  ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J)                        INT02500
  394.       IF (ICC.NE.0) GO TO 1090                                          INT02510
  395.       WRITE(NPRT,9104)                                                  INT02520
  396. 9104  FORMAT(38H0RESULTS OF TRIANGLE TERM IN DIST FUNC)                 INT02530
  397.       GO TO 1100                                                        INT02540
  398. 1090  IF(ICC.GT.0) GO TO 1112                                           INT02550
  399.       WRITE(NPRT,9105)                                                  INT02560
  400. 9105  FORMAT(1H0,'RESULTS USING PRECEDING TAU AT PRESENT CONC')         INT02570
  401.       GO TO 1100                                                        INT02580
  402. 1112  WRITE(NPRT,9102)                                                  INT02590
  403. 9102  FORMAT(1H0,'RESULTS USING PREVIOUS TAU STORED ON DISK')           INT02600
  404. 1100  WRITE(NPRT,9106)                                                  INT02610
  405. 9106  FORMAT(6X,4HDELA,9X,4HDELB,9X,4HDIFF,10X,3HOSM,7X,7HXENERGY,5X,   INT02620
  406.      *7HXVOLUME,4X,4HITER,6X,4HBETA,7X,6HCONVGA,7X,6HCONVGB)            INT02630
  407.       GO TO 3010                                                        INT02640
  408.  
  409.      
  410. C***************** HNC ***********************************************
  411. C     SOLVE HNC INTEGRAL EQUATION BY ITERATION                         *INT02660
  412. 2010  WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG                         INT02680
  413.       WRITE(NPRT,5555)ZE                                                INT02690
  414.       IF(ZE.EQ.0.0) GO TO 2020                                          INT02700
  415. 5555  FORMAT(1X,F10.5)                                                  INT02710
  416.       WRITE(NPRT,9108)                                                  INT02720
  417.       WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS)    INT02730
  418.       WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4                              INT02740
  419. 2020  WRITE(NPRT,9111)                                                  INT02750
  420. 9111  FORMAT(19H0RESULTS OF HNC EQN)                                    INT02760
  421.       WRITE(NPRT,9106)                                                  INT02770
  422. 2030  CONTINUE                                                          INT02780
  423.       DO 2040 I = 1,NIONS                                               INT02790
  424.       DO 2040 J = 1,NIONS                                               INT02800
  425. 2040  CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1)                         INT02810
  426.       DO 2050 I = 1,NIONS                                               INT02820
  427.       DO 2050 J = 1,NIONS                                               INT02830
  428. 2050  CALL FFT3S(F0,ZETA(1,I,J),M1,RL,TZETA(1,I,J),1)                   INT02840
  429.       DO 2070 I = 1,NIONS                                               INT02850
  430.       DO 2070 J = 1,NIONS                                               INT02860
  431.       DO 2070 K = 1,N                                                   INT02870
  432.       SUM1 = 0.0                                                        INT02880
  433.       DO 2060 L = 1,NIONS                                               INT02890
  434. 2060  SUM1 = TZETA(K,I,L)*TH(K,L,J)*C(L)+SUM1                           INT02900
  435. 2070  P1(K,I,J) = SUM1                                                  INT02910
  436.       DO 2090 I = 1,NIONS                                               INT02920
  437.       DO 2090 J = 1,NIONS                                               INT02930
  438.       DO 2090 K = 1,N                                                   INT02940
  439.       P2 = 0.0                                                          INT02950
  440.       P3 = 0.0                                                          INT02960
  441.       DO 2080 L = 1,NIONS                                               INT02970
  442.       P2 = TQ(K,I,L)*TZETA(K,L,J)*C(L)+P2                               INT02980
  443. 2080  P3 = TQ(K,I,L)*P1(K,L,J)*C(L)+P3                                  INT02990
  444. 2090  TTAU(K,I,J) = P1(K,I,J)+P2+P3                                     INT03000
  445.       DO 2100 I = 1,NIONS                                               INT03010
  446.       DO 2100 J = 1,NIONS                                               INT03020
  447. 2100  CALL FFT3S(F0,TTAU(1,I,J),M1,RL,TAU(1,I,J),-1)                    INT03030
  448.       I1 = I1+1                                                         INT03040
  449. 2110  DO 2120 I = 1,NIONS                                               INT03050
  450.       DO 2120 J = 1,NIONS                                               INT03060
  451. C     WRITE(7,5556)                                                     INT03070
  452. 5556  FORMAT(1X,'HI')                                                   INT03080
  453.       DO 2115 K = 1,N                                                   INT03090
  454. C somethign about this raises the floating point exception
  455. C looking deeper, the math library doesnt like that something here has
  456. C an inf value. if is a valid float though. line 39 of
  457. C math_errf.c library.
  458. C it probably would preffer IMPLICIT NONE
  459. C      
  460. 2115  H(K,I,J) = Q1(K,I,J)*EXP(TAU(K,I,J))-1.                           INT03100
  461.       H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5                       INT03110
  462.       DO 2120 K = 1,N                                                   INT03120
  463. 2120  ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J)                        INT03130
  464.       IF (I1.EQ.MAX) GO TO 3010                                         INT03140
  465.       IF ((I1/I2)*I2.EQ.I1) GO TO 3010                                  INT03150
  466.       GO TO 2030
  467.      
  468. C**************** XEnergy*********************************************
  469. C     DELA,DELB,OSM,XENERGY,DFX/DP                                     *INT03180
  470.  
  471. 3010  CONTAC = 0.0                                                      INT03200
  472.       SQ1 = 0.0                                                         INT03210
  473.       SQ2 = 0.0                                                         INT03220
  474.       REMAIN = 0.0                                                      INT03230
  475.       XE = 0.0                                                          INT03240
  476.       DO 3020 I = 1,NIONS                                               INT03250
  477.       DO 3020 J = 1,NIONS                                               INT03260
  478.       NA = MINT(I,J)                                                    INT03270
  479.       HNA = H(NA,I,J)                                                   INT03280
  480.       H(NA,I,J) = HNA*2.+1.                                             INT03290
  481.       AA = A(I,J)**3                                                    INT03300
  482.       NB = NA +NINT(I,J)                                                INT03310
  483.       BB = (A(I,J)+BA(I,J))**3                                          INT03320
  484.       NC = NB+LINT(I,J)                                                 INT03330
  485.       AC = (A(I,J)+BA(I,J)+CB(I,J))**3                                  INT03340
  486.       AD = PI4*(A1(I,J)-A(I,J))                                         INT03350
  487.       CALL TRAP(NA,NB,I,J,SUM1,SUM2,SUM3)                               INT03360
  488.       T(I,J) = PI4*SUM1*SPAC1                                           INT03370
  489.       G(I,J) = PI4*SUM2*SPAC1                                           INT03380
  490.       G4(I,J) = PI4*SUM3*SPAC1                                          INT03390
  491.       GM(I,J) = G(I,J)                                                  INT03400
  492.       TM(I,J) = T(I,J)                                                  INT03410
  493.       XE = XE + FNU(I)*FNU(J)*E(I,J)*(PI43*(BB-AA)+G(I,J))              INT03420
  494.       CALL TRAP(NB,NC,I,J,SUM1,SUM2,SUM3)                               INT03430
  495.       T(I,J) = PI4*SUM1*SPAC1                                           INT03440
  496.       G(I,J) = PI4*SUM2*SPAC1                                           INT03450
  497.       G4(I,J) = PI4*SUM3*SPAC1+G4(I,J)                                  INT03460
  498.       GW(I,J) = G(I,J)                                                  INT03470
  499.       TW(I,J) = T(I,J)                                                  INT03480
  500.       XE = XE+FNU(I)*FNU(J)*W(I,J)*(PI43*(AC-BB)+G(I,J))                INT03490
  501.       ERROR1 = -Z(I)*Z(J)*FY+AD*H(NA,I,J)*XR(NA)                        INT03500
  502.       ERROR2 = -Z(I)*Z(J)*FX+AD*H(NA,I,J)*XR(NA)**2                     INT03510
  503.       ERROR3 = -Z(I)*Z(J)*FZ+AD*H(NA,I,J)*XR(NA)**4                     INT03520
  504. C      PRINT *, "*****From xEnergy*******"
  505. C      PRINT *, "NC",NC,"N",N,"I",I,"J",J
  506. C      PRINT *,"SUM1",SUM1,"SUM2",SUM2,"SUM3",SUM3
  507.  
  508.       CALL SIMP(NC,N,I,J,SUM1,SUM2,SUM3)                                INT03530
  509.       T(I,J) = PI43*SPAC1*SUM1-PI2*A(I,J)**2+TW(I,J)+TM(I,J)+ERROR1     INT03540
  510.       G(I,J) = PI43*(SPAC1*SUM2-A(I,J)**3)+GM(I,J)+GW(I,J)+ERROR2       INT03560
  511.       G4(I,J) = PI43*SPAC1*SUM3-2.513234*A(I,J)**5+G4(I,J)+ERROR3       INT03570
  512.       GA(I,J) = H(NA,I,J)+1.                                            INT03580
  513.       GB(I,J) = H(NB,I,J)+1.                                            INT03590
  514.       GC(I,J) = H(NC,I,J)+1.                                            INT03600
  515.       CONTAC = CONTAC + FNU(I)*FNU(J)*GA(I,J)*AA                        INT03610
  516.       SQ1 = SQ1+FNU(I)*FNU(J)*(GB(I,J)*BB*(1.-EXP(-E(I,J)+W(I,J))))     INT03620
  517.       SQ2 = SQ2+FNU(I)*FNU(J)*(GC(I,J)*AC*(1.-EXP(-W(I,J))))            INT03630
  518.       REMAIN = REMAIN + Z(I)*Z(J)*FNU(I)*FNU(J)*T(I,J)                  INT03640
  519.       IF(I1.EQ.MAX) GO TO 3011                                          INT03650
  520.       H(NA,I,J) = HNA                                                   INT03660
  521.       GO TO 3020                                                        INT03670
  522.  3011 DC(NA,I,J) = DC(NA,I,J)+HNA+1.                                    INT03680
  523.  3020 CONTINUE                                                          INT03690
  524.       CONTAC = CONTAC*2.09439*CC/FNO**2                                 INT03700
  525.       SQ1 = SQ1*2.09439*CC/FNO**2                                       INT03710
  526.       SQ2 = SQ2*2.09439*CC/FNO**2                                       INT03720
  527.       REMAIN = REMAIN*FL1*CC/(6.*FNO**2)                                INT03730
  528.       OSM = 1.+CONTAC+SQ1+SQ2+REMAIN                                    INT03740
  529.       XE = 1.97827*TEMP*(3.*(1.+DT)*REMAIN*FNO+XE*CC*0.5/FNO)           INT03750
  530.       XFP = -8.31467*1.E-05*TEMP*DP*3.*REMAIN*FNO                       INT03760
  531.       GM12 = (G(1,2)+G(2,1))*.5                                         INT03770
  532.       CONVGB = 100*(G(1,2)-G(2,1))/GM12                                 INT03780
  533.       IF (ZE.EQ.0.0) GO TO 3050                                         INT03790
  534.       S2 = 0.0                                                          INT03800
  535.       S4 = 0.0                                                          INT03810
  536.       DO 3040 I = 1,NIONS                                               INT03820
  537.       SUM1 = 0.0                                                        INT03830
  538.       SUM2 = 0.0                                                        INT03840
  539.       DO 3030 J = 1,NIONS                                               INT03850
  540.       SUM1 = SUM1 + C(J)*Z(J)*G(I,J)                                    INT03860
  541. 3030  SUM2 = SUM2 + C(J)*Z(J)*G4(I,J)                                   INT03870
  542.       F2(I) = SUM1                                                      INT03880
  543.       F4(I) = SUM2                                                      INT03890
  544.       Y(I) = 6.*Z(I)/FKFK                                               INT03900
  545.       R2(I) = 1.+F2(I)/Z(I)                                             INT03910
  546.       R4(I) = 1. + F4(I)/Y(I)                                           INT03920
  547.       S2 = C(I)*Z(I)*F2(I)+S2                                           INT03930
  548. 3040  S4 = C(I)*Z(I)*F4(I)+S4                                           INT03940
  549.       CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12 - C(1)+C(2)               INT03950
  550.       CS2 = FKFK/FL                                                     INT03960
  551.       CS4 = 6./FL                                                       INT03970
  552.       E2 = 1.+S2/CS2                                                    INT03980
  553.       E4 = 1.+S4/CS4                                                    INT03990
  554.       BETA = FNO*R2(1)*R2(2)/(FNU(1)*R2(2)+FNU(2)*R2(1))                INT04000
  555. 3050  CONTINUE                                                          INT04010
  556.       FG = 1.0/(CC*GM12)                                                INT04020
  557.       DELA = (FG-1.0-FG*BETA)/(1.0+FG*BETA)                             INT04030
  558.       DELB = FG-1.0                                                     INT04040
  559.       DIFF = DELA-DELB                                                  INT04050
  560.       ITER = I1                                                         INT04060
  561.       FITER = ITER                                                      INT04070
  562.       WRITE(NPRT,9112) DELA,DELB,DIFF,OSM,XE,XFP,ITER,BETA,CONVGA,CONVGBINT04080
  563. 9112  FORMAT(  6E13.4,I5,3E13.4)                                        INT04090
  564.       IF (I1.EQ.MAX) GO TO 3060                                         INT04100
  565.       IF(I1.EQ.0) GO TO 2010                                            INT04110
  566.       IF ((I1/I2)*I2.EQ.I1) GO TO 2030                                  INT04120
  567. 3060  WRITE(NPRT,9113)                                                  INT04130
  568. 9113  FORMAT(15H0MATRIX  G(I,J),21X,15HMATRIX  GM(I,J),20X,             INT04140
  569.      *15HMATRIX  GW(I,J),20X,15HMATRIX  G4(I,J))                        INT04150
  570.       WRITE(NPRT,9114) ((G(I,J),J=1,NIONS),(GM(I,L),L=1,NIONS),(GW(I,M),INT04160
  571.      *M=1,NIONS),(G4(I,K),K=1,NIONS),I=1,NIONS)                         INT04170
  572. 9114  FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6,6X,2E14.6)                   INT04180
  573.       WRITE(NPRT,9115)                                                  INT04190
  574. 9115  FORMAT(15H0MATRIX  T(I,J),21X,15HMATRIX  TM(I,J),20X,             INT04200
  575.      *15HMATRIX  TW(I,J))                                               INT04210
  576.       WRITE(NPRT,9116)((T(I,J),J=1,NIONS),(TM(I,L),L=1,NIONS),(TW(I,M),MINT04220
  577.      *=1,NIONS),I=1,NIONS)                                              INT04230
  578. 9116  FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6)                             INT04240
  579.       WRITE(NPRT,9117)                                                  INT04250
  580. 9117  FORMAT(16H0MATRIX  GA(I,J),20X,15HMATRIX  GB(I,J),20X,            INT04260
  581.      *15HMATRIX  GC(I,J))                                               INT04270
  582.       WRITE(NPRT,9118) ((GA(I,J),J=1,NIONS),(GB(I,L),L=1,NIONS),(GC(I,M)INT04280
  583.      *,M=1,NIONS),I=1,NIONS)                                            INT04290
  584. 9118  FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6)                             INT04300
  585.       WRITE(NPRT,9119)                                                  INT04310
  586. 9119  FORMAT(17H0MATRIX MINT(I,J),3X,16HMATRIX NINT(I,J),3X,            INT04320
  587.      *16HMATRIX LINT(I,J))                                              INT04330
  588.       WRITE(NPRT,9120) ((MINT(I,J),J=1,NIONS),(NINT(I,L),L=1,NIONS),(LININT04340
  589.      *T(I,M),M=1,NIONS),I=1,NIONS)                                      INT04350
  590. 9120  FORMAT(1X,2I4,13X,2I4,13X,2I4)                                    INT04360
  591.       WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG                         INT04370
  592. 9107  FORMAT(8H0CONTAC=,1P,E12.4,3X,'SQUARE1=',1P,E12.4,3X,'SQUARE2=',  INT04380
  593.      *1P,E12.4,3X,'REMAIN=',1P,E12.4,3X,'GAMMA=',1P,E12.4)              INT04390
  594.       IF (ZE.EQ.0.0) GO TO 3070                                         INT04400
  595.       WRITE(NPRT,9108)                                                  INT04410
  596. 9108  FORMAT(6X,1HZ,13X,2HF2,12X,2HR2,12X,1HY,13X,2HF4,12X,2HR4)        INT04420
  597.       WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS)    INT04430
  598. 9109  FORMAT(1X,6E14.6)                                                 INT04440
  599.       WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4                              INT04450
  600. 9110  FORMAT('S2=',E14.6,2X,'CS2=',1P,E14.6,2X,'E2=',1P,E14.6,2X,'S4='  INT04460
  601.      *,1P,E14.6,2X,'CS4=',1P,E14.6,2X,'E4=',1P,E14.6)                   INT04470
  602. 3070  CONTINUE                                                          INT04480
  603.       WRITE(NPRT,9121)                                                  INT04490
  604. 9121  FORMAT(7X,5HXR(N),7X,9H Q(1,1,N),5X,8HH(1,1,N),5X,9H Q(1,2,N),    INT04500
  605.      *6X,8HH(1,2,N),5X,9H Q(2,2,N),6X,8HH(2,2,N))                       INT04510
  606.       WRITE(NPRT,9122) XR(N),Q(N,1,1),H(N,1,1),Q(N,1,2),H(N,1,2),       INT04520
  607.      *Q(N,2,2),H(N,2,2)                                                 INT04530
  608. 9122  FORMAT(7E14.4//)                                                  INT04540
  609.       IF(EX.NE.2.0.AND.EX.NE.3.0) GO TO 4070                            INT04550
  610. 3080  WRITE(NPRT,9123)                                                  INT04560
  611. 9123  FORMAT(7X,5HXR(K),6X,10H DC(1,1,K),5X,8HH(1,1,K),5X,10H DC(1,2,K),INT04570
  612.      *5X,8HH(1,2,K),5X,10H DC(2,1,K),5X,8HH(2,1,K),5X,10H DC(2,2,K),    INT04580
  613.      *5X,8HH(2,2,K))                                                    INT04590
  614.       WRITE(NPRT,9124)(XR(K),((DC(K,I,J),H(K,I,J),J=1,NIONS),I=1,NIONS) INT04600
  615.      *,K=1,N)                                                           INT04610
  616.       WRITE(NDSK,9125)C1,C2                                             INT04620
  617.       WRITE(NDSK,9133)FO,N,M1,RL,((MINT(I,J),I=1,2),J=1,2)              INT04630
  618. 9133  FORMAT(F8.4,2I4,F8.4,4I15.6)                                      INT04640
  619.       WRITE(NDSK,9126)(XR(K),H(K,1,1),H(K,1,2),K=1,N)                   INT04650
  620.       WRITE(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=1,N)       INT04660
  621.       REWIND NDSK2                                                      INT04670
  622. 9124  FORMAT(1X,1P,9E14.4)                                              INT04680
  623. C    ERIC some issue with this line format
  624. 9125  FORMAT(1X,A8,2E16.6)                                              INT04690
  625. 9126  FORMAT(2E16.7)                                                    INT04700
  626.       IF (ZE.EQ.0.0) GO TO 4030                                         INT04710
  627.       IF (EX.NE.3.0) GO TO 4070                                          INT04720
  628.  
  629. C **********************************************************************INT04730
  630. C     CALC OF CHARGE DENSITY                                           *INT04740
  631. C **********************************************************************INT04750
  632.       DO 4020 K = 1,N                                                   INT04760
  633.       DO 4020 I = 1,NIONS                                               INT04770
  634.       SUM1 = 0.0                                                        INT04780
  635.       SUM2 = 0.0                                                        INT04790
  636.       DO 4010 J = 1,NIONS                                               INT04800
  637.       SUM1 = SUM1 + FNU(J)*Z(J)*H(K,I,J)                                INT04810
  638. 4010  SUM2 = SUM2 + FNU(J)*Z(J)*Q(K,I,J)                                INT04820
  639.       RHO(I,K) = SUM1*XR(K)**2                                          INT04830
  640. 4020  RQO(I,K) = SUM2*XR(K)**2                                          INT04840
  641.       GO TO 4060                                                        INT04850
  642. 4030  Z(1) = 1.0                                                        INT04860
  643.       Z(2) = -1.0                                                       INT04870
  644.       DO 4050 K = 1,N                                                   INT04880
  645.       DO 4050 I = 1,NIONS                                               INT04890
  646.       SUM1 = 0.0                                                        INT04900
  647.       DO 4040 J = 1,NIONS                                               INT04910
  648. 4040  SUM1 = SUM1+FNU(J)*Z(J)*H(K,I,J)                                  INT04920
  649.       RHO(I,K) = SUM1*XR(K)**2                                          INT04930
  650. 4050  RQO(I,K) = 0.0                                                    INT04940
  651. 4060  WRITE(NPRT,9127)                                                  INT04950
  652. 9127  FORMAT('1CHARGE DENSITY'///1X,9X,'XR',11X,'RHO1',10X,'RQO1',11X,  INT04960
  653.      *'RHO2',10X,'RQO2')                                                INT04970
  654. C      PRINT *, 'NCD:' , NCD
  655. C      PRINT *, 'K:' , K
  656. C      PRINT *, 'IC:' , IC
  657.       WRITE(NPRT,9128)(XR(K),RHO(1,K),RQO(1,K),RHO(2,K),RQO(2,K),K=1,NCDINT04980
  658.      *,IC)                                                              INT04990
  659.       PRINT *, 'NR:' , NR
  660.  
  661.  
  662. 9128  FORMAT(1X,1P,5E14.4)                                              INT05000
  663.       WRITE(NDSK,9126)(RHO(1,K),RHO(2,K),K=1,N2)                        INT05010
  664. 9129  FORMAT(2E16.7)                                                    INT05020
  665.       REWIND NDSK                                                       INT05030
  666. 4070  CONTINUE                                                          INT05040
  667.       B(1) = 0.0                                                        INT05050
  668.       B(2) = 1.0                                                        INT05060
  669.       Z(1) = ZA(1)                                                      INT05070
  670.       Z(2) = ZA(2)                                                      INT05080
  671.       ZE = Z(1)                                                         INT05090
  672.       IF (ZE.EQ.0.0) GO TO 5040                                         INT05100
  673. C **********************************************************************INT05110
  674. C     CALC OF DIELECTRIC RESPONSE FUNCTION                              INT05120
  675. C **********************************************************************INT05130
  676.       DO 5010 I = 1,NIONS                                               INT05140
  677.       DO 5010 J = 1,NIONS                                               INT05150
  678.       H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5                       INT05160
  679. 5010  CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1)                         INT05170
  680.       SPACT = PI2/(2.*FNN-1.)/SPAC1                                     INT05180
  681.       TK = 1.0                                                          INT05190
  682.       DO 5030 K = 1,N                                                   INT05200
  683.       Y1(K) = TK*SPACT                                                  INT05210
  684.       SUM = 0.0                                                         INT05220
  685.       DO 5020 I = 1,NIONS                                               INT05230
  686.       DO 5020 J = 1,NIONS                                               INT05240
  687. 5020  SUM = SUM + C(I)*C(J)*Z(I)*Z(J)*TH(K,I,J)                         INT05250
  688.       Y1K2 = Y1(K)**2                                                   INT05260
  689.       DRF(K) = 1.-(FKFK/Y1K2)*(1.+SUM/FC)                               INT05270
  690.       DRDF(K) = 1.-FKFK/(FKFK+Y1K2)                                     INT05280
  691. 5030  TK = TK+1.                                                        INT05290
  692.       WRITE(NPRT,9130)                                                  INT05300
  693. 9130  FORMAT('1DIELECTRIC RESPONSE FUNCTION'///1X,3(9X,'YK',12X,'DRF',  INT05310
  694.      *11X,'DRDF'))                                                      INT05320
  695.       WRITE(NPRT,9131)(Y1(K),DRF(K),DRDF(K),K=1,NR,IR)                  INT05330
  696.  
  697.      
  698. 9131  FORMAT(1X,1P,9E14.4)                                              INT05340
  699. 5040  CONTINUE
  700. C5040 CALL TIMER(DATE,TIME,IVIRT,ITOT)                                  INT05350
  701. C      RVIRT=FLOAT(IVIRT)*1.E-6                                          INT05360
  702. C      WRITE(NPRT,9132)TIME,DATE,RVIRT                                   INT05370
  703.       WRITE(NPRT,9132)date_string                                        INT05370
  704. C9132  FORMAT(//////////'END OF EXECUTION AT',A8,' ON ',A8,5X,F8.3,      INT05380
  705. C     *' SECONDS USED')                                                  INT05390
  706. 9132  FORMAT(//////////'END EXECUTION ON ',A20)                                                  INT05390
  707.  
  708.       IF (ABS(CHARGE-0.0).LT.1.E-04) GO TO 120                          INT05400
  709. C      GOTO 40                                                           INT05410
  710. C     Program loop      
  711. 5050  CONTINUE                                                          INT05420
  712.       STOP                                                              INT05430
  713.       END                                                               INT05440
  714. C ****Program stops and ends ****************************************
  715.  
  716. C************************ SUBROUTINES *******************************
  717.  
  718. C************************ ADJ ***************************************
  719.       SUBROUTINE ADJ(A,A1,M,SPAC1)                                      INT05450
  720.       S1 = A/SPAC1                                                      INT05460
  721.       L = S1+1                                                          INT05470
  722.       S2 = L                                                            INT05480
  723.       IF(S2-S1.GE.0.5) GO TO 10                                         INT05490
  724.       A1 = S2*SPAC1                                                     INT05500
  725.       M = L                                                             INT05510
  726.       GO TO 20                                                          INT05520
  727. 10    A1 = (S2-1.)*SPAC1                                                INT05530
  728.       M = L-1                                                           INT05540
  729. 20    CONTINUE                                                          INT05550
  730.       RETURN                                                            INT05560
  731.       END
  732.      
  733. C************************ TRAP ***************************************
  734.       SUBROUTINE TRAP(N1,N11,I,J,SUM1,SUM2,SUM3)                        INT05580
  735.       DIMENSION H(2048,2,2),XR(2048)                                    INT05590
  736.       COMMON/JAY/ H,XR                                                  INT05600
  737.       IF (N11.EQ.N1) GO TO 20                                           INT05610
  738.       SUM1 = (H(N1,I,J)*XR(N1)+H(N11,I,J)*XR(N11))*.5                   INT05620
  739.       SUM2 = (H(N1,I,J)*XR(N1)**2+H(N11,I,J)*XR(N11)**2)*.5             INT05630
  740.       SUM3 = (H(N1,I,J)*XR(N1)**4+H(N11,I,J)*XR(N11)**4)*.5             INT05640
  741.       L1 = N1+1                                                         INT05650
  742.       L2 = N11-1                                                        INT05660
  743.       IF (N11.EQ.L1) GO TO 30                                           INT05670
  744.       DO 10 K = L1,L2                                                   INT05680
  745.       SUM1 = SUM1 + H(K,I,J)*XR(K)                                      INT05690
  746.       SUM3 = SUM3 + H(K,I,J)*XR(K)**4                                   INT05700
  747. 10    SUM2 = SUM2 + H(K,I,J)*XR(K)**2                                   INT05710
  748.       GO TO 30                                                          INT05720
  749. 20    SUM1 = 0.0                                                        INT05730
  750.       SUM2 = 0.0                                                        INT05740
  751.       SUM3 = 0.0                                                        INT05750
  752. 30    RETURN                                                            INT05760
  753.       END
  754. C********************** SIMPSON **************************************
  755.       SUBROUTINE SIMP(NC,N,I,J,SUM1,SUM2,SUM3)                          INT05780
  756.       DIMENSION H(2048,2,2),XR(2048)                                    INT05790
  757.       COMMON /JAY/ H,XR                                                 INT05800
  758.  
  759.       SUM1 = H(NC,I,J)*XR(NC)                                           INT05810
  760.       SUM2 = SUM1*XR(NC)                                                INT05820
  761.       SUM3 = SUM2*XR(NC)**2                                             INT05830
  762. C      PRINT *, 'NC:',NC
  763. C     Setting NC to 1 for testing??
  764.       NC = 1
  765.       N2 = (N-NC-2)/2        
  766. C      PRINT *, '**N2:', N2      
  767.       DO 10 K = 1,N2                                                    INT05850
  768.       K4 = NC+2*K-1                                                     INT05860
  769.       K5 = K4+1                                                         INT05870
  770.       K6 = NC+2*N2+1                                                    INT05880
  771.       K7 = K6+1                                                         INT05890
  772.       SUM1 = SUM1 + 4.*H(K4,I,J)*XR(K4)+2.*H(K5,I,J)*XR(K5)             INT05900
  773.       SUM3 = SUM3 + 4.*H(K4,I,J)*XR(K4)**4+2.*H(K5,I,J)*XR(K5)**4       INT05910
  774. 10    SUM2 = SUM2 + 4.*H(K4,I,J)*XR(K4)**2+2.*H(K5,I,J)*XR(K5)**2       INT05920
  775.       SUM1 = SUM1 + 4.*H(K6,I,J)*XR(K6)+2.*H(K7,I,J)*XR(K7)             INT05930
  776.       SUM2 = SUM2 + 4.*H(K6,I,J)*XR(K6)**2+2.*H(K7,I,J)*XR(K7)**2       INT05940
  777.       SUM3 = SUM3 + 4.*H(K6,I,J)*XR(K6)**4+2.*H(K7,I,J)*XR(K7)**4       INT05950
  778.       RETURN                                                            INT05960
  779.       END
  780. C**************** Fast Fourier Transform *************************      
  781.       SUBROUTINE FFT3S(F0,F,MM,RL,T,MODE)                               INT05990
  782. C     FAST FOURIER TRANSFORM OF 3-DIMENSIONAL CENTROSYMMETRIC FUNCTION  INT06000
  783.       DIMENSION F(2048),T(1)                                               INT06010
  784.       COMMON /FDATA/ N,N2,FF,RTTWO,NC,NS,ND,NR                          INT06020
  785.       COMMON /FWORK/ W(12288)                                           INT06030
  786. C     DIMENSION OF W SHOULD BE .GE. 3*N                                 INT06040
  787. C     FIXED DATA PREPARATION                                            INT06050
  788.       DATA PI/3.141593/,M/-1/                                           INT06060
  789.       MD = MM+1                                                         INT06070
  790.       CALL SETFT(MD,M)                                                  INT06080
  791.       N = 2**MM                                                         INT06090
  792.       RN = FLOAT(N)                                                     INT06100
  793.       DR = RL/RN                                                        INT06110
  794.       RL3 = RL*RL*RL                                                    INT06120
  795.       C = 4.*RL3*FF/RN                                                  INT06130
  796.       CC = 0.0                                                          INT06140
  797.       A = 1.0                                                           INT06150
  798.       IF(IABS(MODE)-2)202,200,201                                       INT06160
  799. 200   CC = 0.5*PI*DR*DR*DR*F(1)                                         INT06170
  800.       A = .125*RN*RN*RN/(RL3*RL3)                                       INT06180
  801.       GOTO 202                                                          INT06190
  802. 201   CC = F0                                                           INT06200
  803. 202   IF(MODE.GT.0)GOTO 203                                             INT06210
  804.       C = 1.0/C                                                         INT06220
  805.       CC = A*CC                                                         INT06230
  806. 203   CONTINUE                                                          INT06240
  807. C     COPY INPUT VECTOR INTO WORKING STORAGE                            INT06250
  808.       DO 204 I=1,N                                                      INT06260
  809. 204   W(ND+I+1) = C*FF*F(I)*I                                           INT06270
  810. C     Eric Starts debugging here
  811. C     RETURN
  812. C     END
  813.       CALL RSA(M)                                                       INT06280
  814. C     COPY INTO OUTPUT VECTOR
  815. C     N eq 2048      
  816. C      DO 205 I=1,N    
  817.       DO 205 I=1,1                                                  
  818. 205   T(I) = W(NR+I+1)/FLOAT(I)+CC                                      INT06310
  819.       RETURN                                                            INT06320
  820.       END                                                               INT06330
  821. C**************************** SETFT **********************************
  822.       SUBROUTINE SETFT(MM,M)                                            INT06340
  823. C     PREPARE FIXED DATA FOR FAST FOURIER TRANSFORM                     INT06350
  824.       COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES                      INT06360
  825.       COMMON /FWORK/ W(12288)                                           INT06370
  826.       IF(MM.LT.0)RETURN                                                 INT06380
  827.       IF(MM.EQ.M)RETURN                                                 INT06390
  828.       M = MM                                                            INT06400
  829.       N = 2**M                                                          INT06410
  830.       N2 = N/2                                                          INT06420
  831.       F = 1.0/FLOAT(N)                                                  INT06430
  832.       DA = 6.2831853071796D0*F                                          INT06440
  833.       RTTWO = SQRT(2.0)                                                 INT06450
  834.       F = SQRT(F)                                                       INT06460
  835.       N1 = (N2-1)/2                                                     INT06470
  836.       NC = 0                                                            INT06480
  837.       NS = N1                                                           INT06490
  838.       NDATA = N1+N1                                                     INT06500
  839.       IF(N1.LE.0)RETURN                                                 INT06510
  840.       DO 10 I=1,N1                                                      INT06520
  841.       A = FLOAT(I)*DA                                                   INT06530
  842.       W(I) = COS(A)                                                     INT06540
  843. 10    W(N1+I) = SIN(A)                                                  INT06550
  844.       RETURN                                                            INT06560
  845.       END  
  846. C************************** RSA *************************************  
  847.       SUBROUTINE RSA(M)                                                 INT06580
  848. C     REAL ODD ANALYSIS OR SYNTHESIS                                    INT06590
  849.       COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES                      INT06600
  850.       COMMON /FWORK/ W(12288)                                           INT06610
  851.       NQ = N2                                                           INT06620
  852.       NIB = NDATA +1                                                    INT06630
  853.       NOB = NDATA+N2+2                                                  INT06640
  854.       DO 200 L = 1,M                                                    INT06650
  855.       NIM = NIB+NQ                                                      INT06660
  856.       NIE = NIM+NQ                                                      INT06670
  857.       NQH = NQ/2                                                        INT06680
  858.       NQL = (NQ-1)/2                                                    INT06690
  859.       NOM = NOB+NQL+1                                                   INT06700
  860. C     TPRIME = 0                                                        INT06710
  861.       IF(NQL.EQ.0)GOTO 20                                               INT06720
  862.       DO 10 IT = 1,NQL                                                  INT06730
  863.       W(NOB+IT)=W(NIB+IT)-W(NIM-IT)                                     INT06740
  864. 10    W(NOM+IT) = W(NIB+IT)+W(NIM-IT)                                   INT06750
  865. 20    IF(NQH.NE.NQL) W(NOM) = 2.0*W(NIB+NQH)
  866. C ??  0.LT.TPRIME.LT.NP/2.                                              INT06770
  867.       NO1 = NQ                                                          INT06780
  868.       NO2 = N2-NQ                                                       INT06790
  869.       IF(NO1-NO2)40,80,120                                              INT06800
  870. 40    NI=NO1+NO1                                                        INT06810
  871.       W(NOB+NO1) = W(NIB+NI)+W(NIM+NI)                                  INT06820
  872.       W(NOB+NO2) = -W(NIB+NI)+W(NIM+NI)                                 INT06830
  873.       IF(NQL.EQ.0)GOTO 60                                               INT06840
  874.       CC = W(NC+NO1)                                                    INT06850
  875.       SS = W(NS+NO1)                                                    INT06860
  876.       DO 50 IT = 1,NQL                                                  INT06870
  877.       RE = CC*W(NIE+NI-IT)-SS*W(NIM+NI-IT)                              INT06880
  878.       AI = SS*W(NIE+NI-IT)+CC*W(NIM+NI-IT)                              INT06890
  879.       W(NOM+NO1+IT) = W(NIM+NI+IT)-RE                                   INT06900
  880.       W(NOM+NO2+IT) = W(NIM+NI+IT)+RE                                   INT06910
  881.       W(NOB+NO1+IT) = W(NIB+NI+IT)+AI                                   INT06920
  882. 50    W(NOB+NO2+IT) = -W(NIB+NI+IT)+AI                                  INT06930
  883. 60    IF(NQH.EQ.NQL)GOTO 70                                             INT06940
  884.       NR = NO1/2                                                        INT06950
  885.       W(NOM+NO1) = 2.0*(W(NS+NR)*W(NIM+NI+NQH)+W(NC+NR)*W(NIB+NI+NQH))  INT06960
  886.       W(NOM+NO2) = 2.0*(W(NC+NR)*W(NIM+NI+NQH)-W(NS+NR)*W(NIB+NI+NQH))  INT06970
  887. 70    NO1 = NO1+NQ                                                      INT06980
  888.       NO2 = NO2-NQ                                                      INT06990
  889.       IF(NO1-NO2)40,80,120                                              INT07000
  890. 80    W(NOB+NO1)=W(NIM)                                                 INT07010
  891.       IF(NQL.EQ.0)GOTO 100                                              INT07020
  892.       DO 90 IT = 1,NQL                                                  INT07030
  893.       W(NOM+NO1+IT) = W(NIM+IT)                                         INT07040
  894. 90    W(NOB+NO1+IT) = W(NIE-IT)                                         INT07050
  895. 100   IF(NQH.NE.NQL)W(NOM+NO1) = RTTWO*W(NIM+NQH)                       INT07060
  896. 120   NT = NIB                                                          INT07070
  897.       NIB = NOB                                                         INT07080
  898.       NOB = NT                                                          INT07090
  899.       NQ = NQH                                                          INT07100
  900. 200   CONTINUE                                                          INT07110
  901.       NRES = NIB-1                                                      INT07120
  902.       RETURN                                                            INT07130
  903.       END                                                               INT07140
  904.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!