PROGRAM GAMMA1 !*********************************************************************************** ! This program uses the sigma profiles of two pure components to calculate ! the liquid-phase activity coefficients in a solution. This is the first ! step in predicting VLE for mixtures. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! THIS PROGRAM WRITTEN BY: ! RICHARD OLDLAND (roldland@vt.edu) MIKE ZWOLAK (zwolak@caltech.edu) ! DEPARTMENT OF CHEMICAL ENGINEERING PHYSICS DEPARTMENT ! VIRGINIA TECH CALIFORNIA INSTITUTE OF TECHNOLOGY ! BLACKSBURG, VA 24060 PASADENA, CA 91125 ! ! PHYSICAL CONSTANTS AND PARAMETERS: ! EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom) ! AEFFPRIME = EFFECTIVE SURFACE AREA (ANGSTROMS**2) --FROM LIN ! RGAS = IDEAL GAS CONSTANT (Kcal/mol*K) ! VNORM = VOLUME NORMALIZATION CONSTANT (A**3) --FROM LIN ! ANORM = AREA NORMALIZATION CONSTANT (A**2) --FROM LIN ! COORD = THE COORIDINATION NUMBER --FROM LIN ! CHB = HYDROGEN BONDING COEFFICIENT (Kcal/mole*Angstroms**4/e**2) ! SIGMAHB = CUTOFF VALUE FOR HYDROGEN BONDING (e/Angstrom**2) ! EPS = RELATIVE PERMITTIVITY --FROM LIN ! ALPHAPRIME = A CONSTANT USED IN THE MISFIT ENERGY CALCULATION ! ! INPUT PARAMETERS: ! SYSTEMP = THE SYSTEM TEMPERATURE (K) ! COMP = NUMBER OF COMPONENTS IN THE SYSTEM --SET TO 2 FOR BINARY ! SYSCOMP = NAMES OF COMPONENTS IN THE SYSTEM ! VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3) ! ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM ! OF THE INDIVIDUAL PROFILE. ! ! ! LITERATURE CITED: ! Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the ! Quantitative Calculation of Solvation Phenomena. J. Phys. Chem 1995, 99, 2224. ! Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. Refinement and Parameterization of ! COSMO-RS. J. Phys. Chem A 1998, 102, 5074. ! Klamt, A.; Eckert, F.; COSMO-RS: A Novel and Efficient Method for the a Priori ! Prediction of Thermophysical Data of Liquids. Fluid Phase Equilibria 2000, ! 172, 43. ! Lin, S.T.; Sandler, S. A Priori Phase Equilibrium Prediction from a Segment ! Contribution Solvation Model. Ind. Eng. Chem. Res, 2002, 41, 899 ! Lin, S.T.; Quantum Mechanical Approaches to the Prediction of Phase Equilibria: ! Solvation Thermodynamics and Group Contribution Methods, PhD. Dissertation, ! University of Delaware, Newark, DE, 2000 ! ! ! PROGRAM CURRENTLY SETUP FOR BINARY MIXTURES ONLY !********************************************************************************** IMPLICIT NONE REAL,PARAMETER:: EO = 2.395*10.0**-4, AEFFPRIME = 7.5, RGAS = 0.001987 REAL,PARAMETER:: VNORM = 66.69, ANORM = 79.53 REAL :: FPOL, ALPHA, ALPHAPRIME, COORD, EPS, SYSTEMP, SIGMAHB, CHB, FRAC1, FRAC2 REAL :: SYSPRES, SIGMAACC, SIGMADON, SUMMATION, BOTTHETA, BOTPHI, PHI1, PHI2 REAL :: THETA1, THETA2, L1, L2, GAMMASG1, GAMMASG2, GAMMA, GAMMA2, SUMGAMMA1 REAL :: SUMGAMMA2, N1, N2, LNGAMMA, LNGAMMA2 INTEGER :: I, J, K, L, M, COMPSEG, COUNT, COMP REAL, DIMENSION(2):: VCOSMO, ACOSMO, RNORM, QNORM REAL, DIMENSION(:), ALLOCATABLE:: COUNTER, DENOM, PROFILE, NUMER, SEGGAMMA REAL, DIMENSION(:), ALLOCATABLE:: SEGGAMMAOLD, CONVERG REAL, DIMENSION(:,:), ALLOCATABLE:: SIGMA, DELTAW, SEGGAMMAPR, SEGGAMMAOLDPR, CONPR CHARACTER (16), DIMENSION(2):: SYSCOMP CHARACTER (256), DIMENSION(2):: FILENAME CHARACTER :: ANSWER CHARACTER (256):: OUTPUT COMPSEG = 51 !NUMBER OF INTERVALS FOR THE SIGMA PROFILE EPS = 3.667 !(LIN AND SANDLER USE A CONSTANT FPOL WHICH YEILDS EPS=3.68) COORD = 10.0 !(KLAMT USED 7.2) SIGMAHB = 0.0084 CHB = 85580.0 COMP =2 FPOL = (EPS-1.0)/(EPS+0.5) ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ALPHAPRIME = FPOL*ALPHA ALLOCATE(SIGMA(COMPSEG,COMP), COUNTER(COMPSEG), PROFILE(COMPSEG), NUMER(COMPSEG),& DENOM(COMPSEG), DELTAW(COMPSEG,COMPSEG), SEGGAMMA(COMPSEG), SEGGAMMAOLD(COMPSEG), & CONVERG(COMPSEG), SEGGAMMAPR(COMPSEG,COMP), SEGGAMMAOLDPR(COMPSEG,COMP), & CONPR(COMPSEG,COMP)) !DEFINE SYSTEM TEMPERATURE (K) WRITE(*,*) "ENTER IN THE SYSTEM TEMPERATURE (K)" READ(*,*) SYSTEMP !DEFINE THE SYSTEM AS WELL AS THE AREA AND VOL FROM THE COSMO CALCULATION DO I=1, COMP WRITE(*,*) "ENTER THE NAME OF COMPONENT", I READ(*,*) SYSCOMP(I) WRITE(*,*) "ENTER THE LOCATION AND NAME OF THE ", SYSCOMP(I), "SIGMA PROFILE (256 character max)" READ(*,*) FILENAME(I) WRITE(*,*) "ENTER THE CAVITY VOLUME FOR COMPONENT",I READ(*,*) VCOSMO(I) END DO !OPEN THE SIGMA PROFILES FOR THE PURE COMPONENTS OPEN(UNIT=1, FILE=FILENAME(1), STATUS="OLD", ACTION="READ", POSITION="REWIND") OPEN(UNIT=2, FILE=FILENAME(2), STATUS="OLD", ACTION="READ", POSITION="REWIND") !READ INDIVIDUAL SIGMA PROFILES; COUNTER IS THE SIGMA VALUE, SIGMA IS P(SIGMA) DO K = 1, COMP ACOSMO(K) =0.0 DO J=1, COMPSEG READ(K,*)COUNTER(J), SIGMA(J,K) ACOSMO(K)=ACOSMO(K)+SIGMA(J,K) END DO END DO CLOSE(1) CLOSE(2) !ESTABLISH OUTPUT FILES. THE PROGRAM CREATES A FILE FOR THE GAMMA-X DATA WRITE(*,*) "ENTER THE LOCATION AND NAME FOR THE OUTPUT FILE, INCLUDE EXTENSION." READ (*,*) OUTPUT OPEN(UNIT=11, FILE = OUTPUT, STATUS="NEW") WRITE(11,*) "SYSTEMP", SYSTEMP, "KELVIN" 5 FORMAT (1X,A10,5X,A12,5X,A12,5X,A12,5X,A12) 6 FORMAT (1X,A10,5X,A9,5X,A10,7X,A10,5X,A10) WRITE(11,6) "MOLE FRAC", "GAMMA1", "GAMMA2", "LNGAMMA1", "LNGAMMA2" WRITE(11,5) "X1", SYSCOMP(1), SYSCOMP(2), SYSCOMP(1), SYSCOMP(2) !VARYING MOLE FRACTIONS // ONLY WORKS FOR BINARY MIXTURE DO FRAC1 = 0.005, 0.995, 0.01 FRAC2 = 1.0 - FRAC1 !CALCULATE THE MIXTURE SIGMA PROFILE DO J =1,COMPSEG NUMER(J) = FRAC1*SIGMA(J,1) + FRAC2*SIGMA(J,2) DENOM(J) = FRAC1*ACOSMO(1) + FRAC2*ACOSMO(2) PROFILE(J) = NUMER(J)/DENOM(J) END DO DO I = 1, COMPSEG DO K = 1, COMPSEG IF (COUNTER(I)>=COUNTER(K)) THEN SIGMAACC = COUNTER(I) SIGMADON = COUNTER(K) END IF IF (COUNTER(I)