DIFFERENT PROGRAMMES COSMO6SAC ET NRTL-SAC I/COSMO-SAC I.1/ORIGINAL COSMO-SAC (2002) I.1.1/COSMO-SAC binaire FORTRAN I.1.2/COSMO-SAC binaire MATLAB I.1.3/COSMO-SAC ternaire FORTRAN I.2/COSMO-SAC (2007) I.2.1/COSMO-SAC binaire I.2.1.1/COSMO-SAC binaire FORTRAN I.2.1.2/COSMO-SAC binaire MATLAB I.2.2/COSMO-SAC ternaire I.3/COSMO-SAC (2010) I.3.1/COSMO-SAC binaire I.3.1.1/COSMO-SAC binaire FORTRAN I.3.1.2/COSMO-SAC binaire MATLAB I.3.2/COSMO-SAC ternaire II/NRTL-SAC II.1/NRTL-SAC binaire II.1.1/NRTL-SAC binaire FORTRAN II.1.2/NRTL-SAC binaire MATLAB II.2/NRTL-SAC ternaire I/COSMO-SAC I.1/ORIGINAL COSMO-SAC (2002) I.1.1/COSMO-SAC binaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PROGRAM WRITTEN BY: ! ERIC MULLINS (PMULLINS@VT.EDU) ! DEPARTMENT OF CHEMICAL ENGINEERING ! VIRGINIA TECH ! BLACKSBURG, VA 24061 ! ! MODIFIED BY BAPTISTE BOUILLOT ! LABORATOIRE DE GENIE CHIMIQUE ! TOULOUSE FRANCE ! ! 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. ! ! IN DATA.txt FILE !!!!!!!!! ! ! 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 PURE SOLVENT SOLUBILITY PREDICTIONS ONLY ! TO PREDICT BINARY SOLVENT SOLUBILITY, MODIFY THE MIXTURE SIGMA PROFILE (LINE 245) ! MODIFY READ INPUT FILE (LINE 130) ! MODIFY NUMBER OF COMPONENTS, VARIABLE: [COMP], CURRENTLY SET TO 2 !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.5D0 ! ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA, DELTAW, WHB, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE, SIGMA1SOLVENT, COUNTERSOLUTE !SIGMA2SOLVENT, REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONPR REAL*8, DIMENSION (:), ALLOCATABLE :: SEGGAMMA, SEGGAMMAOLD, CONVERG, X1DATA, X2DATA, X !, X3DATA, REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, COUNTER, MIXPROFILE, DELTAHF, Cp REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, SIGMAHB, CHB, SIGMAACC, SIGMADON, SIGMAHBNEW, SUMMATION REAL*8 :: RELTOL, CHANGE, X2OLD, X1NEW, X2NEW, MFSUM ! X3NEW INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMTEMP INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU ! UNIT2SOLV, CHARACTER (128) :: INPUTFILE, OUTPUTFILE CHARACTER (16), DIMENSION (:), ALLOCATABLE :: SYSCOMP CHARACTER (50), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE !, SOLVENT2FILE, CHARACTER (16), DIMENSION (:), ALLOCATABLE :: SOLV1 !, SOLV2 CHARACTER (15), DIMENSION (:), ALLOCATABLE :: SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 SIGMAHB = 0.0084D0 ! HYDROGEN-BONDING INTERACTION CUT-OFF, e/ANG^2 SIGMAHBNEW = 0.022D0 ! FOR USE WITH NEW EXPRESSION FOR EXCHANGE ENERGY, e/ANG^2 CHB = 85580.0D0 ! HYDROGEN-BONDING INTERACTION CONSTANT, KCAL*ANG^4/MOL/e^2 FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 2 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data-Cp.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMTEMP = 41 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMTEMP, NUMPOINTS CLOSE(12) write(*,*) NUMTEMP, NUMPOINTS ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), & VCOSMO2SOLVENT(NUMPOINTS), & SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS), TM(NUMPOINTS), DELTAHF(NUMPOINTS),& SYSTEMP(NUMPOINTS), Cp(NUMPOINTS), & X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS),SOLUTEFILE(NUMPOINTS), SOLVENT1FILE(NUMPOINTS), & UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS), SIGMASOLUTE(COMPSEG,NUMPOINTS), & COUNTERSOLUTE(COMPSEG,NUMPOINTS), ACOSMOSOLUTE(NUMPOINTS), & SIGMA1SOLVENT(COMPSEG,NUMPOINTS), COUNTER1SOLVENT(COMPSEG,NUMPOINTS), STAT = ISTAT) !, & !SIGMA2SOLVENT(COMPSEG,NUMPOINTS), COUNTER2SOLVENT(COMPSEG,NUMPOINTS), & !ACOSMO2SOLVENT(NUMPOINTS)), SOLV2(NUMPOINTS), !UNIT2SOLV(NUMPOINTS),X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), SOLVENT2FILE(NUMPOINTS)) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMTEMP, NUMPOINTS DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), TM(L), DELTAHF(L), & SYSTEMP(L), Cp(L), X2DATA(L) ! , X1DATA(L), X3DATA(L), SOLV2(L), VCOSMO2SOLVENT(L) END DO CLOSE (12) !################################### ! LA MIENNE : ! ! DO L=1, NUMPOINTS ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") ! READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), TM(L), DELTAHF(L), & ! SYSTEMP(L), X2DATA(L) ! , X1DATA(L), X3DATA(L), SOLV2(L), VCOSMO2SOLVENT(L) ! CLOSE(12) ! END DO ! #################### REPRISE : ################### DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNITSOLU (I) = I + 10020 SOLUTEFILE (I) = ""//SOLU(I) ![USER SPECIFIED FILE LOCATION] SOLVENT1FILE (I) = ""//SOLV1(I) ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLUTEFILE(I), X2DATA(I), TM(I) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT(J,K) !write(*,*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT(J,K), ACOSMO1SOLVENT(K) !test test test END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A7,1X,A7) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","GAMMA1","GAMMA2" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP), COUNTER(COMPSEG), MIXPROFILE(COMPSEG), DELTAW(COMPSEG,COMPSEG),& SEGGAMMA(COMPSEG),SEGGAMMAOLD(COMPSEG),CONVERG(COMPSEG),SEGGAMMAPURE(COMPSEG,COMP),& SEGGAMMAPUREOLD(COMPSEG,COMP), CONVPURE(COMPSEG,COMP), WHB(COMPSEG, COMPSEG), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS !write(*,*) SYSTEMP(P) ! BEGIN ITERATION LOOP FOR ALL POINTS (température), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMTEMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 WHB = 0.0D0 SIGMAACC = 0.0D0 SIGMADON = 0.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 DO I = 1, COMPSEG SIGMA(I,1) = SIGMA1SOLVENT(I,P) SIGMA(I,2) = SIGMASOLUTE(I,P) COUNTER(I) = COUNTER1SOLVENT(I,P) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMOSOLUTE(P) X(1) = 1.0D0 - X2DATA(P) X(2) = X2DATA(P) !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) X2OLD = X(2) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] MIXPROFILE(J) = ((X(1)*SIGMA(J,1))+(X(2)*SIGMA(J,2))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) ![MIXTURE PROFILE FOR BINARY SOLVENT SYSTEM] !MIXPROFILE(J)=((X(1)*SIGMA(J,1))+(X(2)*SIGMA(J,2))+(X(3)*SIGMA(J,3))) / & !((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))+(X(3)*ACOSMO(3))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO I = 1, COMPSEG DO K = 1, COMPSEG ! ASSIGN HYDROGEN BONDING DONOR AND ACCEPTOR IF (COUNTER(I)>=COUNTER(K)) THEN SIGMAACC = COUNTER(I) SIGMADON = COUNTER(K) END IF IF (COUNTER(I)= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(2) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(2) = (X(2) + 2.0D0*X2OLD)/3.0D0 MFSUM = X(1)+X(2) ! +X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM ! Pour afficher la fraction molaire finale intermédiaire ! write(*,*) 'final', X(2) ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(2) - X2OLD)/(X2OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),GAMMA(1),GAMMA(2) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(2) write (16,*) GAMMA(2) write(*,*) X(2) ! #### CAS OU ON INCREMENTE LA TEMPERATURE : SYSTEMP(P) = SYSTEMP(P) + 1 END DO ! Iteration temperature (variable P) END DO ! Iteration points (variable R) CLOSE (13) ! CLOSE OUTPUT FILE close (15) ! CLOSE l'AUTRE OUTPUT FILE close (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, WHB, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, SYSTEMP, & X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE, COUNTERSOLUTE, ACOSMOSOLUTE, SIGMA1SOLVENT, COUNTER1SOLVENT, & ACOSMO1SOLVENT, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" !SOLV2, VCOSMO2SOLVENT, SOLVENT2FILE, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT !X1DATA, X3DATA, UNIT2SOLV, END PROGRAM GAMMASOLUBILITY I.1.2/COSMO-SAC 2002 MATLAB function [y] = cosmosac02() % % cosmosac02erreur renvoie la solubilité prédite par cosmosac de 2002 avec pour entrée % Tm la température de fusion, Hf, l'enthalpie de fusion, n la taille du vecteur Tm ou Hf, et p le numéro du mélange à traiter %format long; % %i %chb=X(1); %S0=X(2); comp=2; % compseg = 51;% NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) EO = 0.0002395 ;% PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) R = 0.001987 ;% KCAL/(MOL K) aeffprime = 7.50 ;% ANGSTROMS^2 ANORM = 79.5300 ;% NORMALIZED SURFACE AREA, ANGSTROMS^2 VNORM = 66.6900 ;% NORMALIZED CAVITY VOLUME, ANGSTROMS^3 COORD = 10.0000 ;% COORDINATE NUMBER, Z, KLAMT SET TO 7.2 EPS = 3.667 ;% DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 CHB = 85580.0; FPOL = (EPS-1.0)/(EPS+0.5); ALPHA = (0.3*aeffprime^(1.5))/(EO) ; alphaprime = FPOL*ALPHA ; sigmahb = 0.0084 ;% Hydrogen bonding interaction cut-off SIGMAHBNEW = 0.022 ; % % Lecture des données d'entrée et initialisation de la boucle A = importdata('input02.txt'); % A(p,1)=solu,2=solv,3=Vcosmosolu,4=vcosmosolv,5=Tm,6=Hm,7=T,8=X2 % A.textdata{p,1}=solu.txt, (2)=solv.txt %A.data(p,1)=Vcosmosolu, (p,2)=Vcosmosolv, (p,3)=Tm, (p,4)=Hm, (p,5)=T, (p,6)=X2 N = size(A.data,1); sol=zeros(N,1); solub=zeros(N,1); mixprofile=zeros(compseg,2); x=zeros(2,1); sigmasol=zeros(compseg,2); retol = 1E-6; %sigmasolv=zeros(compseg,2,2); % Solute = load(A.textdata{1,1}); Solv = load(A.textdata{1,2}); for p=1:N % ajouter une condition sur le cas précédent pour pas tout refaire... % w=p-1; % if (p~=1&&A.textdata{p,1}~=A.textdata{w,1}) %aime pas le (p-1) en argument... Solute = load(A.textdata{p,1}); % end Solv = load(A.textdata{p,2}); Delta=zeros(51,51); Acosmo=zeros(2); Acosmo(1)=ones(1,51)*Solv(:,2); Acosmo(2)=ones(1,51)*Solute(:,2); x(2)=A.data(p,6); x(1)=1-x(2); iter=-1; %while err > 10^-8 % boucle sur l'erreur entre solubilité prédite et exp change = 1; while change > retol % boucle sur la détermination de x(2) iter = iter +1 ; % Calcul du gamma combinatoire lngammasg = lngammac(A.data(p,1),A.data(p,2),Acosmo(1),Acosmo(2),x(2),VNORM,ANORM,COORD); % % Calcul du gamma résiduel sigmasol(:,1)=Solv(:,2); sigmasol(:,2)=Solute(:,2); mixprofile=((x(1)*sigmasol(:,1))+(x(2)*sigmasol(:,2)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); %x(1)=1-x(2); x2old=x(2); for i=1:compseg for k=1:compseg if (Solv(i,1)>=Solv(k,1)) sigmaacc = Solv(i,1); sigmadon = Solv(k,1); end if (Solv(i,1)=0.00001 % boucle pour seggamma du mélange seggammaold = seggamma; for i=1:compseg summation = 0; for k=1:compseg summation = summation + mixprofile(k)*seggammaold(k)*... exp(-Delta(i,k)/(R*A.data(p,5))); end % en retirant la boucle sur k: summation = summation + mixprofile(:,h)'*(seggammaold(:,h).*exp(-Delta(i,:,b,h)'./(R*A.data(p,5)))); % malheureusement, c'est moins intéressant niveau optimisation. pkoi ??? seggamma(i)=exp(-log(summation)); seggamma(i)=(seggamma(i)+2*seggammaold(i))/3; end converg=abs((seggamma-seggammaold)./seggammaold); end % Fin boucle if iter==0 conpr=ones(51,2); seggammapureold=ones(compseg,comp); for l=1:comp seggammapure=ones(compseg,comp); while (max(max(max(conpr(:,l))))>=0.00001) % début boucle pour calcul de seggammapure conpr(:,l)=zeros(51,1); %zeros(51,2); seggammapureold(:,l)=seggammapure(:,l); for i=1:compseg summation=0; for k=1:51 summation=summation+(sigmasol(k,l)/Acosmo(l))... *seggammapureold(k,l)*exp(-Delta(i,k)/(R*A.data(p,5))); end seggammapure(i,l)=exp(-log(summation)); seggammapure(i,l)=(seggammapure(i,l)+seggammapureold(i,l))/2; end conpr(:,l)=abs((seggammapure(:,l)-seggammapureold(:,l))./seggammapureold(:,l)); end %fin boucle seggammapure (WHILE) end %fin boucle sur l end % % Calcul de gamma residuel gamma=ones(2,1); lngamma=ones(2,1); sumgamma=zeros(2,1); for k=1:2 for i=1:compseg sumgamma(k)=sumgamma(k) + ((sigmasol(i,k)/aeffprime)*(log(seggamma(i)/... seggammapure(i,k)))); end gamma(k)=exp(sumgamma(k)+(lngammasg(k))); lngamma(k)=log(gamma(k)); end % Recalcul de la fraction molaire x(2) = exp(A.data(p,4)/8.3145./A.data(p,3)*(1-A.data(p,3)./A.data(p,5))-lngamma(2)); x(2) = (x(2)+2*x2old)/3; mfsum=x(1)+x(2); x(1)=x(1)/mfsum; x(2)=x(2)/mfsum; change=abs((x(2)-x2old)/(x2old+1E-16)); end disp(lngammasg(2)); disp(sumgamma(2)); sol(p)=x(2); end y=sol; end function [y]=lngammac(Vcosmosolute,Vcosmosolv,Acosmosolv,Acosmosolute,Xs,Vnorm,Anorm,COORD) %theta=zeros(2,1); %phi=zeros(2,1); %lsg=zeros(2,1); %lngammasg=zeros(2,1); X(1)=1-Xs; X(2)=Xs; Rnorm(1)=Vcosmosolv/Vnorm; Rnorm(2)=Vcosmosolute/Vnorm; Qnorm(1)=Acosmosolv/Anorm; Qnorm(2)=Acosmosolute/Anorm; theta=(X.*Qnorm)./(X*Qnorm'); phi=(X.*Rnorm)./(X*Rnorm'); lsg=(COORD/2).*(Rnorm-Qnorm)-(Rnorm-1); lngammasg=log(phi./X)+(COORD/2).*Qnorm.*log(theta./phi)+lsg... - (phi./X)*(X*lsg'); %[X(2),Rnorm(2), Qnorm(2), theta(2), phi(2), lsg(2)] y=lngammasg; end I.1.3/COSMO-SAC ternaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PROGRAM WRITTEN BY: ! ERIC MULLINS (PMULLINS@VT.EDU) ! DEPARTMENT OF CHEMICAL ENGINEERING ! VIRGINIA TECH ! BLACKSBURG, VA 24061 ! ! MODIFIED BY BAPTISTE BOUILLOT ! LABORATOIRE DE GENIE CHIMIQUE ! TOULOUSE FRANCE ! ! 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. ! ! IN DATA.txt FILE !!!!!!!!! ! ! 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 PURE SOLVENT SOLUBILITY PREDICTIONS ONLY ! TO PREDICT BINARY SOLVENT SOLUBILITY, MODIFY THE MIXTURE SIGMA PROFILE (LINE 245) ! MODIFY READ INPUT FILE (LINE 130) ! MODIFY NUMBER OF COMPONENTS, VARIABLE: [COMP], CURRENTLY SET TO 2 !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.5D0 ! ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA, DELTAW, WHB, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE, SIGMA1SOLVENT, COUNTERSOLUTE, SIGMA2SOLVENT REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONPR REAL*8, DIMENSION (:), ALLOCATABLE :: SEGGAMMA, SEGGAMMAOLD, CONVERG, X1DATA, X2DATA, X , X3DATA REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, COUNTER, MIXPROFILE, DELTAHF, DELTACP REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, SIGMAHB, CHB, SIGMAACC, SIGMADON, SIGMAHBNEW, SUMMATION REAL*8 :: RELTOL, CHANGE, X2OLD, X1NEW, X2NEW, MFSUM , X3OLD INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMTEMP INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU, UNIT2SOLV CHARACTER (128) :: INPUTFILE, OUTPUTFILE CHARACTER (16), DIMENSION (:), ALLOCATABLE :: SYSCOMP CHARACTER (50), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE, SOLVENT2FILE CHARACTER (16), DIMENSION (:), ALLOCATABLE :: SOLV1 , SOLV2 CHARACTER (15), DIMENSION (:), ALLOCATABLE :: SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 SIGMAHB = 0.0084D0 ! HYDROGEN-BONDING INTERACTION CUT-OFF, e/ANG^2 SIGMAHBNEW = 0.022D0 ! FOR USE WITH NEW EXPRESSION FOR EXCHANGE ENERGY, e/ANG^2 CHB = 85580.0D0 ! HYDROGEN-BONDING INTERACTION CONSTANT, KCAL*ANG^4/MOL/e^2 FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 3 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data_ternaire_cp.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMTEMP = 41 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMTEMP, NUMPOINTS CLOSE(12) write(*,*) NUMTEMP, NUMPOINTS ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), & VCOSMO2SOLVENT(NUMPOINTS), DELTACP(NUMPOINTS), & SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS), TM(NUMPOINTS), DELTAHF(NUMPOINTS), & SYSTEMP(NUMPOINTS), X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS),SOLUTEFILE(NUMPOINTS), & SOLVENT1FILE(NUMPOINTS), UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS), & SIGMASOLUTE(COMPSEG,NUMPOINTS), COUNTERSOLUTE(COMPSEG,NUMPOINTS), & ACOSMOSOLUTE(NUMPOINTS), SIGMA1SOLVENT(COMPSEG,NUMPOINTS), & COUNTER1SOLVENT(COMPSEG,NUMPOINTS), SIGMA2SOLVENT(COMPSEG,NUMPOINTS), & COUNTER2SOLVENT(COMPSEG,NUMPOINTS), ACOSMO2SOLVENT(NUMPOINTS), SOLV2(NUMPOINTS), & UNIT2SOLV(NUMPOINTS),X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), SOLVENT2FILE(NUMPOINTS), & STAT = ISTAT) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMTEMP, NUMPOINTS DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), SOLV2(L), VCOSMO2SOLVENT(L), & TM(L), DELTAHF(L), DELTACP(L), SYSTEMP(L), X2DATA(L)!, X3DATA(L) END DO CLOSE (12) DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNIT2SOLV (I) = I + 140 UNITSOLU (I) = I + 10020 SOLUTEFILE (I) = ""//TRIM(SOLU(I))//".txt" ![USER SPECIFIED FILE LOCATION] SOLVENT1FILE (I) = ""//TRIM(SOLV1(I))//".txt" ![USER SPECIFIED FILE LOCATION] SOLVENT2FILE (I) = ""//TRIM(SOLV2(I))//".txt" ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLVENT2FILE(I), SOLUTEFILE(I), TM(I) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT(J,K) END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT2SOLV(K), FILE = SOLVENT2FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO2SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT2SOLV(K),*) COUNTER2SOLVENT(J,K), SIGMA2SOLVENT(J,K) ACOSMO2SOLVENT(K) = ACOSMO2SOLVENT(K) + SIGMA2SOLVENT(J,K) END DO CLOSE (UNIT2SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A2) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)", "COMP(3)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","X3" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP), COUNTER(COMPSEG), MIXPROFILE(COMPSEG), DELTAW(COMPSEG,COMPSEG),& SEGGAMMA(COMPSEG),SEGGAMMAOLD(COMPSEG),CONVERG(COMPSEG),SEGGAMMAPURE(COMPSEG,COMP),& SEGGAMMAPUREOLD(COMPSEG,COMP), CONVPURE(COMPSEG,COMP), WHB(COMPSEG, COMPSEG), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS !write(*,*) SYSTEMP(P) ! BEGIN ITERATION LOOP FOR ALL POINTS (température), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMTEMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 WHB = 0.0D0 SIGMAACC = 0.0D0 SIGMADON = 0.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 DO I = 1, COMPSEG SIGMA(I,1) = SIGMA1SOLVENT(I,P) SIGMA(I,2) = SIGMA2SOLVENT(I,P) SIGMA(I,3) = SIGMASOLUTE(I,P) COUNTER(I) = COUNTER1SOLVENT(I,P) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMO2SOLVENT(P) ACOSMO(3) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMO2SOLVENT(P) VCOSMO(3) = VCOSMOSOLUTE(P) X(3) = 0.000000001D0 X(2) = X2DATA(P) X(1) = 1.0D0 - X2DATA(P) - X3DATA(P) !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) - X(3) X2OLD = X(2) X3OLD = X(3) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] MIXPROFILE(J) = ((X(1)*SIGMA(J,1))+(X(2)*SIGMA(J,2))+(X(3)*SIGMA(J,3))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO I = 1, COMPSEG DO K = 1, COMPSEG ! ASSIGN HYDROGEN BONDING DONOR AND ACCEPTOR IF (COUNTER(I)>=COUNTER(K)) THEN SIGMAACC = COUNTER(I) SIGMADON = COUNTER(K) END IF IF (COUNTER(I)= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(3) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(3) = (X(3) + 2.0D0*X3OLD)/3.0D0 MFSUM = X(1)+X(2)+X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM X(3) = X(3)/MFSUM ! Pour afficher la fraction molaire finale intermédiaire ! write(*,*) 'final', X(2) ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(3) - X3OLD)/(X3OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),X(3) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(3) write (16,*) GAMMA(3) write(*,*) X(3) ! #### CAS OU ON INCREMENTE LA TEMPERATURE : SYSTEMP(P) = SYSTEMP(P) + 1 END DO ! Iteration temperature (variable P) END DO ! Iteration points (variable R) CLOSE (13) ! CLOSE OUTPUT FILE close (15) ! CLOSE l'AUTRE OUTPUT FILE close (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, WHB, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, SYSTEMP, DELTACP, & X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE, COUNTERSOLUTE, ACOSMOSOLUTE, SIGMA1SOLVENT, COUNTER1SOLVENT, & ACOSMO1SOLVENT, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT, & SOLV2, UNIT2SOLV,X1DATA, X3DATA, SOLVENT2FILE, & VCOSMO2SOLVENT, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" !SOLV2, VCOSMO2SOLVENT, SOLVENT2FILE, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT !X1DATA, X3DATA, UNIT2SOLV, END PROGRAM GAMMASOLUBILITY I.2/COSMO-SAC (2007) I.2/COSMO-SAC binaire I.2.1/COSMO-SAC binaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PHYSICAL CONSTANTS AND PARAMETERS: ! EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom) ! S0 = SHAPE PARAMETER ? ! 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) ! 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 ! VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3) ! ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM ! OF THE INDIVIDUAL PROFILE. ! ! IN data_refined_Cp.txt FILE !!!!!!!!! ! ! PROGRAM CURRENTLY SETUP FOR PURE SOLVENT SOLUBILITY PREDICTIONS ONLY !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.25D0 ! ou 7.5 selon publi ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 REAL*8, PARAMETER :: S0 = 0.007D0 ! shape parameter INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:,:,:), ALLOCATABLE :: DELTAW REAL*8, DIMENSION (:,:,:), ALLOCATABLE :: SIGMA, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONPR REAL*8, DIMENSION (:,:), ALLOCATABLE :: WHB, SEGGAMMA, SEGGAMMAOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE_NHB, SIGMASOLUTE_HB, MIXPROFILE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_HB, COUNTERSOLUTE REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONVERG REAL*8, DIMENSION (:), ALLOCATABLE :: X1DATA, X2DATA, X , COUNTER!, X3DATA, REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, DELTAHF, DELTACP REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, CHB, SUMMATION REAL*8 :: RELTOL, CHANGE, X2OLD, X1NEW, X2NEW, MFSUM ! X3NEW INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMTEMP, H, B INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU ! UNIT2SOLV, CHARACTER (40) :: INPUTFILE, OUTPUTFILE CHARACTER (30), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE !, SOLVENT2FILE, CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLV1 !, SOLV2 CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 CHB = 3484.42D0 !28476.21D0 !(2004) HYDROGEN-BONDING INTERACTION CONSTANT, KCAL*ANG^4/MOL/e^2 !FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS FPOL = 0.6916D0 !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 2 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data_refined_Cp.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMTEMP = 41 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMTEMP, NUMPOINTS CLOSE(12) write(*,*) NUMTEMP, NUMPOINTS ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), & VCOSMO2SOLVENT(NUMPOINTS), SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS),& TM(NUMPOINTS), DELTAHF(NUMPOINTS), SYSTEMP(NUMPOINTS), DELTACP(NUMPOINTS), & X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS),SOLUTEFILE(NUMPOINTS), SOLVENT1FILE(NUMPOINTS), & UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS), SIGMASOLUTE_NHB(COMPSEG,NUMPOINTS), & SIGMASOLUTE_HB(COMPSEG,NUMPOINTS), COUNTERSOLUTE(COMPSEG,NUMPOINTS), ACOSMOSOLUTE(NUMPOINTS), & SIGMA1SOLVENT_NHB(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_HB(COMPSEG,NUMPOINTS), & COUNTER1SOLVENT(COMPSEG,NUMPOINTS), STAT = ISTAT) !, & !SIGMA2SOLVENT(COMPSEG,NUMPOINTS), COUNTER2SOLVENT(COMPSEG,NUMPOINTS), & !ACOSMO2SOLVENT(NUMPOINTS)), SOLV2(NUMPOINTS), !UNIT2SOLV(NUMPOINTS),X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), SOLVENT2FILE(NUMPOINTS)) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMTEMP, NUMPOINTS DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), TM(L), DELTAHF(L), & SYSTEMP(L), DELTACP(L), X2DATA(L) ! , X1DATA(L), X3DATA(L), SOLV2(L), VCOSMO2SOLVENT(L) END DO CLOSE (12) !################################### ! LA MIENNE : ! ! DO L=1, NUMPOINTS ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") ! READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), TM(L), DELTAHF(L), & ! SYSTEMP(L), X2DATA(L) ! , X1DATA(L), X3DATA(L), SOLV2(L), VCOSMO2SOLVENT(L) ! CLOSE(12) ! END DO ! #################### REPRISE : ################### DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNITSOLU (I) = I + 10020 SOLVENT1FILE (I) = TRIM(SOLV1(I))//'_refined.txt' ![USER SPECIFIED FILE LOCATION] SOLUTEFILE (I) = TRIM(SOLU(I))//'_refined.txt' ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLUTEFILE(I), X2DATA(I), TM(I) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT_NHB(J,K), SIGMA1SOLVENT_HB(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT_NHB(J,K) + SIGMA1SOLVENT_HB(J,K) !write(*,*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT(J,K), ACOSMO1SOLVENT(K) !test test test END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE_NHB(J,K), SIGMASOLUTE_HB(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE_NHB(J,K) + SIGMASOLUTE_HB(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins_refined.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_refined_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_refined_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A7,1X,A7) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","GAMMA1","GAMMA2" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP, 2), COUNTER(COMPSEG),& MIXPROFILE(COMPSEG,2), DELTAW(COMPSEG,COMPSEG,2,2),& SEGGAMMA(COMPSEG,2),SEGGAMMAOLD(COMPSEG,2),CONVERG(COMPSEG,2),SEGGAMMAPURE(COMPSEG,COMP,2),& SEGGAMMAPUREOLD(COMPSEG,COMP,2), CONVPURE(COMPSEG,COMP), WHB(COMPSEG, COMPSEG), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP, 2), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS !write(*,*) SYSTEMP(P) ! BEGIN ITERATION LOOP FOR ALL POINTS (température), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMTEMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 WHB = 0.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 DO I = 1, COMPSEG COUNTER(I) = COUNTER1SOLVENT(I,P) !COUNTER(I) = COUNTER1SOLVENT(I,1) ! 1->NHB, 2->HB !COUNTER(I) = COUNTER1SOLVENT(I,2) ! SIGMA(segment, species, nhb or hb) SIGMA(I,1,1) = SIGMA1SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMA1SOLVENT_NHB(I,P) SIGMA(I,2,1) = SIGMASOLUTE_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMASOLUTE_NHB(I,P) SIGMA(I,1,2) = SIGMA1SOLVENT_HB(I,P)*(1-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,2,2) = SIGMASOLUTE_HB(I,P)*(1-dexp(-COUNTER(I)**2/(2*S0**2))) !write(*,*) (dexp(-COUNTER(I)**2/(2*S0**2))) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMOSOLUTE(P) X(1) = 1.0D0 - X2DATA(P) X(2) = X2DATA(P) !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) X2OLD = X(2) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] ! MIXPROFILE IS p(sigmai(J)) FOR THE MIXTURE in case NHB (1), HB (2) MIXPROFILE(J,1) = ((X(1)*SIGMA(J,1,1))+(X(2)*SIGMA(J,2,1))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) MIXPROFILE(J,2) = ((X(1)*SIGMA(J,1,2))+(X(2)*SIGMA(J,2,2))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO H=1,2 DO B=1,2 DO I = 1, COMPSEG DO K = 1, COMPSEG IF (H==B.AND.H==2.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN ! DELTAW(I,K,H,B) is delta for segment I for H=1 (NHB) or H=2 (HB) and segment K DELTAW(I,K,H,B) = (ALPHAPRIME/2.0D0)*(COUNTER(I)+COUNTER(K))**2.0D0 - & CHB*(COUNTER(I)-COUNTER(K))**2 ELSE DELTAW(I,K,H,B) = (ALPHAPRIME/2.0D0)*(COUNTER(I)+COUNTER(K))**2.0D0 ENDIF END DO END DO ENDDO ENDDO ! ITERATION FOR MIXTURE SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 SEGGAMMA = 1.0D0 !INTIAL VALUE FOR MIXTURE SEGMENT ACTIVITY COEF. DO SEGGAMMAOLD = SEGGAMMA DO B=1,2 DO I = 1, COMPSEG SUMMATION = 0.0D0 DO H=1,2 DO K = 1, COMPSEG SUMMATION = SUMMATION + MIXPROFILE(K,H)*SEGGAMMAOLD(K,H) * & DEXP(-DELTAW(I,K,B,H)/(RGAS * SYSTEMP(P))) ENDDO ENDDO SEGGAMMA(I,B) = DEXP(-DLOG(SUMMATION)) SEGGAMMA(I,B) = (SEGGAMMA(I,B) + SEGGAMMAOLD(I,B))/2.0D0 END DO ENDDO DO B=1, 2 DO I=1, COMPSEG CONVERG(I,B) = DABS((SEGGAMMA(I,B) - SEGGAMMAOLD(I,B))/SEGGAMMAOLD(I,B)) END DO ENDDO IF (MAXVAL(CONVERG) <=0.000001) EXIT ENDDO ! ITERATION FOR PURE SPECIES, L, SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 ! ITERATION FOR SEGMENT ACITIVITY COEF (PURE SPECIES) ! HERE p(sigma) is just SIGMA/Ai total DO L = 1, COMP SEGGAMMAPURE (:,L,:) = 1.0D0 DO SEGGAMMAPUREOLD (:,L,:) = SEGGAMMAPURE (:,L,:) DO I = 1, COMPSEG DO B=1,2 SUMMATION = 0.0D0 DO H=1,2 DO K = 1, COMPSEG SUMMATION = SUMMATION + (SIGMA(K,L,H)/ACOSMO(L)) * & SEGGAMMAPUREOLD(K,L,H)*DEXP(-DELTAW(I,K,B,H)/(RGAS*SYSTEMP(P))) ENDDO END DO SEGGAMMAPURE(I,L,B)=DEXP(-DLOG(SUMMATION)) SEGGAMMAPURE(I,L,B)=(SEGGAMMAPURE(I,L,B)+SEGGAMMAPUREOLD(I,L,B))/2.0D0 ENDDO END DO DO B=1,2 DO I=1,COMPSEG CONPR(I,L,B)=DABS((SEGGAMMAPURE(I,L,B)-SEGGAMMAPUREOLD(I,L,B))/& SEGGAMMAPUREOLD(I,L,B)) ENDDO ENDDO IF (MAXVAL(CONPR) <=0.000001) EXIT ENDDO ENDDO ! THE STAVERMAN-GUGGENHEIM EQUATION DO I = 1,COMP RNORM(I) = VCOSMO(I)/VNORM QNORM(I) = ACOSMO(I)/ANORM END DO DO I = 1, COMP THETA(I) = (X(I)*QNORM(I))/(X(1)*QNORM(1) + X(2)*QNORM(2)) PHI(I) = (X(I)*RNORM(I))/(X(1)*RNORM(1) + X(2)*RNORM(2)) LSG(I) = (COORD/2.0D0)*(RNORM(I)-QNORM(I))-(RNORM(I)-1.0D0) END DO ! GAMMASG1 AND GAMMASG2 ARE ACTUALLY LNGAMMASG DO I = 1, COMP LNGAMMASG(I) = DLOG(PHI(I)/X(I)) + (COORD/2.0D0)*QNORM(I)* & DLOG(THETA(I)/PHI(I)) + LSG(I) - (PHI(I)/X(I))*(X(1)*LSG(1) + X(2)*LSG(2)) END DO !CALCULATION OF GAMMAS SUMGAMMA = 0.0D0 DO K = 1, COMP DO B = 1, 2 DO I=1, COMPSEG SUMGAMMA(K) = SUMGAMMA(K)+((SIGMA(I,K,B)/AEFFPRIME)*(DLOG(SEGGAMMA(I,B)/ & SEGGAMMAPURE(I,K,B)))) END DO ENDDO GAMMA(K) = DEXP(SUMGAMMA(K) + (LNGAMMASG(K))) LNGAMMA(K) = DLOG(GAMMA(K)) END DO ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(2) = DEXP((DELTAHF(P)/8.3145D0/TM(P))*(1.0D0-TM(P)/SYSTEMP(P))- & (DELTACP(P)/8.3145D0)*(DLOG(TM(P)/SYSTEMP(P))-(TM(P)/SYSTEMP(P))+1.0D0) & - LNGAMMA(2)) ! Affichage des résultats intermédiaires ! WRITE(*,*) P, X(2) !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(2) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(2) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(2) = (X(2) + 2.0D0*X2OLD)/3.0D0 MFSUM = X(1)+X(2) ! +X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM ! Pour afficher la fraction molaire finale intermédiaire ! write(*,*) 'final', X(2) ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(2) - X2OLD)/(X2OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),GAMMA(1),GAMMA(2) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(2) write (16,*) GAMMA(2) write(*,*) X(2) ! #### CAS OU ON INCREMENTE LA TEMPERATURE : SYSTEMP(P) = SYSTEMP(P) + 5 END DO ! Iteration temperature (variable P) END DO ! Iteration points (variable R) CLOSE (13) ! CLOSE OUTPUT FILE close (15) ! CLOSE l'AUTRE OUTPUT FILE close (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, DELTACP, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, WHB, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, SYSTEMP, & X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE_NHB, SIGMASOLUTE_HB, COUNTERSOLUTE, ACOSMOSOLUTE,& SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_HB,COUNTER1SOLVENT, & ACOSMO1SOLVENT, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" !SOLV2, VCOSMO2SOLVENT, SOLVENT2FILE, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT !X1DATA, X3DATA, UNIT2SOLV, END PROGRAM GAMMASOLUBILITY I.2.1.2/COSMO-SAC 2007 binaire MATLAB function [y] = cosmosac(X) %[chb S0]) % % Cosmo-sac renvoie l'erreur entre la solubilité attendue et celle trouvée % on rentre comme parametre : CHB et S0 %format long; % % chb=X(1); S0=X(2); comp=2; compseg = 51;% NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) EO = 0.0002395 ;% PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) R = 0.001987 ;% KCAL/(MOL K) aeffprime = 7.2500 ;% ANGSTROMS^2 ANORM = 79.5300 ;% NORMALIZED SURFACE AREA, ANGSTROMS^2 VNORM = 66.6900 ;% NORMALIZED CAVITY VOLUME, ANGSTROMS^3 COORD = 10.0000 ;% COORDINATE NUMBER, Z, KLAMT SET TO 7.2 EPS = 3.667 ;% DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 FPOL = 0.6916; ALPHA = (0.3000*aeffprime^(1.5))/(EO) ; alphaprime = FPOL*ALPHA ; %SIGMAHB = 0.022 ;% Hydrogen bonding interaction cut-off %SIGMANEW = 0 ; %CHB = 85580 ; %FPOL = (EPS-1.0000)/(EPS+0.5000) ; % % Lecture des données d'entrée et initialisation de la boucle A = importdata('input2.txt'); % A(p,1)=solu,2=solv,3=Vcosmosolu,4=vcosmosolv,5=Tm,6=Hm,7=T,8=X2 % A.textdata{p,1}=solu.txt, (2)=solv.txt %A.data(p,1)=Vcosmosolu, (p,2)=Vcosmosolv, (p,3)=Tm, (p,4)=Hm, (p,5)=T, (p,6)=X2 N = size(A.data,1) sol=zeros(N,1); solub=zeros(N,1); mixprofile=zeros(compseg,2); x=zeros(2,1); sigmasol=zeros(compseg,2,2); retol = 1E-6; %sigmasolv=zeros(compseg,2,2); % Solute = load(A.textdata{1,1}); Solv = load(A.textdata{1,2}); for p=1:N % ajouter une condition sur le cas précédent pour pas tout refaire... % w=p-1; % if (p~=1&&A.textdata{p,1}~=A.textdata{w,1}) %aime pas le (p-1) en argument... Solute = load(A.textdata{p,1}); % end Solv = load(A.textdata{p,2}); Delta=zeros(51,51,2,2); Acosmo=zeros(2); Acosmo(1)=ones(1,51)*Solv(:,2)+ones(1,51)*Solv(:,3); Acosmo(2)=ones(1,51)*Solute(:,2)+ones(1,51)*Solute(:,3); x(2)=A.data(p,6); x(1)=1-x(2); iter=-1; %while err > 10^-8 % boucle sur l'erreur entre solubilité prédite et exp change = 1; while change > retol % boucle sur la détermination de x(2) iter = iter +1 ; % Calcul du gamma combinatoire lngammasg = lngammac(A.data(p,1),A.data(p,2),Acosmo(1),Acosmo(2),x(2),VNORM,ANORM,COORD); % % Calcul du gamma résiduel sigmasol(:,1,1)=Solv(:,2)+Solv(:,3).*(exp(-(Solv(:,1).^2)./(2*S0^2))); sigmasol(:,2,1)=Solute(:,2)+Solute(:,3).*(exp(-(Solute(:,1).^2)./(2*S0^2))); sigmasol(:,1,2)=Solv(:,3).*(1-exp(-(Solv(:,1).^2)./(2*S0^2))); sigmasol(:,2,2)=Solute(:,3).*(1-exp(-(Solute(:,1).^2)./(2*S0^2))); mixprofile(:,1)=((x(1)*sigmasol(:,1,1))+(x(2)*sigmasol(:,2,1)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); mixprofile(:,2)=((x(1)*sigmasol(:,1,2))+(x(2)*sigmasol(:,2,2)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); %x(1)=1-x(2); x2old=x(2); for h=1:2 for b=1:2 for i=1:compseg for k=1:compseg if (h==b&&h==2&&(Solv(i,1)*Solv(k,1))<0) Delta(i,k,h,b)=(alphaprime/2)*(Solv(i,1)+Solv(k,1))^2-chb*(Solv(i,1)-Solv(k,1))^2; else Delta(i,k,h,b)=(alphaprime/2)*(Solv(i,1)+Solv(k,1))^2; end end end end end seggamma = ones(compseg,2); converg=ones(compseg,2); while max(max(converg))>=0.00001 % boucle pour seggamma du mélange seggammaold = seggamma; for b=1:2 for i=1:compseg summation = 0; for h=1:2 for k=1:compseg summation = summation + mixprofile(k,h)*seggammaold(k,h)*... exp(-Delta(i,k,b,h)/(R*A.data(p,5))); end % en retirant la boucle sur k: summation = summation + mixprofile(:,h)'*(seggammaold(:,h).*exp(-Delta(i,:,b,h)'./(R*A.data(p,5)))); % malheureusement, c'est moins intéressant niveau optimisation. pkoi ??? end seggamma(i,b)=exp(-log(summation)); seggamma(i,b)=(seggamma(i,b)+2*seggammaold(i,b))/3; end end converg=abs((seggamma-seggammaold)./seggammaold); end % Fin boucle if iter==0 conpr=ones(51,2,2); seggammapureold=ones(compseg,comp,2); for l=1:comp seggammapure=ones(compseg,comp,2); while (max(max(max(conpr(:,l,:))))>=0.00001) % début boucle pour calcul de seggammapure conpr(:,l,:)=zeros(51,2); %zeros(51,2,2); seggammapureold(:,l,:)=seggammapure(:,l,:); for i=1:compseg for b=1:2 summation=0; for h=1:2 for k=1:51 summation=summation+(sigmasol(k,l,h)/Acosmo(l))... *seggammapureold(k,l,h)*exp(-Delta(i,k,b,h)/(R*A.data(p,5))); end end seggammapure(i,l,b)=exp(-log(summation)); seggammapure(i,l,b)=(seggammapure(i,l,b)+seggammapureold(i,l,b))/2; end end conpr(:,l,:)=abs((seggammapure(:,l,:)-seggammapureold(:,l,:))./seggammapureold(:,l,:)); end %fin boucle seggammapure (WHILE) end %fin boucle sur l end % % Calcul de gamma residuel gamma=ones(2,1); lngamma=ones(2,1); sumgamma=zeros(2,1); for k=1:2 for b=1:2 for i=1:compseg sumgamma(k)=sumgamma(k) + ((sigmasol(i,k,b)/aeffprime)*(log(seggamma(i,b)/... seggammapure(i,k,b)))); end end gamma(k)=exp(sumgamma(k)+(lngammasg(k))); lngamma(k)=log(gamma(k)); end % Recalcul de la fraction molaire x(2) = exp((A.data(p,4)/8.3145/A.data(p,3))*(1-A.data(p,3)/A.data(p,5))-lngamma(2)); x(2) = (x(2)+2*x2old)/3; mfsum=x(1)+x(2); x(1)=x(1)/mfsum; x(2)=x(2)/mfsum; change=abs((x(2)-x2old)/(x2old+1E-16)); end sol(p)=x(2); solub(p)=(log(A.data(p,6))-log(sol(p)))/sqrt(N); end %y=sol; sol y=solub; end function [y]=lngammac(Vcosmosolute,Vcosmosolv,Acosmosolv,Acosmosolute,Xs,Vnorm,Anorm,COORD) %theta=zeros(2,1); %phi=zeros(2,1); %lsg=zeros(2,1); %lngammasg=zeros(2,1); X(1)=1-Xs; X(2)=Xs; Rnorm(1)=Vcosmosolv/Vnorm; Rnorm(2)=Vcosmosolute/Vnorm; Qnorm(1)=Acosmosolv/Anorm; Qnorm(2)=Acosmosolute/Anorm; theta=(X.*Qnorm)./(X*Qnorm'); phi=(X.*Rnorm)./(X*Rnorm'); lsg=(COORD/2).*(Rnorm-Qnorm)-(Rnorm-1); lngammasg=log(phi./X)+(COORD/2).*Qnorm.*log(theta./phi)+lsg... - (phi./X)*(X*lsg'); %[X(2),Rnorm(2), Qnorm(2), theta(2), phi(2), lsg(2)] y=lngammasg; end I.2.2/COSMO-SAC ternaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PHYSICAL CONSTANTS AND PARAMETERS: ! EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom) ! S0 = SHAPE PARAMETER ? ! 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) ! 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 ! VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3) ! ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM ! OF THE INDIVIDUAL PROFILE. ! ! IN DATA.txt FILE !!!!!!!!! ! ! PROGRAM CURRENTLY SETUP FOR PURE SOLVENT SOLUBILITY PREDICTIONS ONLY !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.25D0 ! ou 7.5 selon publi ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 !REAL*8, PARAMETER :: S0 = 0.05D0 !0.007D0 ! shape parameter REAL*8 :: S0,S1 INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:,:,:), ALLOCATABLE :: DELTAW REAL*8, DIMENSION (:,:,:), ALLOCATABLE :: SIGMA, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONPR REAL*8, DIMENSION (:,:), ALLOCATABLE :: SEGGAMMA, SEGGAMMAOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE_NHB, SIGMASOLUTE_HB, MIXPROFILE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_HB, COUNTERSOLUTE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA2SOLVENT_NHB, SIGMA2SOLVENT_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONVERG REAL*8, DIMENSION (:), ALLOCATABLE :: X1DATA, X2DATA, X , COUNTER, X3DATA REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, DELTAHF, DELTACP REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, CHB, SUMMATION REAL*8 :: RELTOL, CHANGE, X1OLD, X2OLD, X3OLD, X1NEW, X2NEW, X3NEW, MFSUM INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMTEMP, H, B INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU , UNIT2SOLV CHARACTER (40) :: INPUTFILE, OUTPUTFILE CHARACTER (30), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE , SOLVENT2FILE CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLV1, SOLV2 CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 !CHB = 3484.42D0 !28476.21D0 !(2004) HYDROGEN-BONDING INTERACTION CONSTANT, KCAL*ANG^4/MOL/e^2 !FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS FPOL = 0.6916D0 !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 3 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data_refined_ternaire_cp.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMTEMP = 41 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMTEMP, NUMPOINTS, CHB, S0, S1 CLOSE(12) write(*,*) NUMTEMP, NUMPOINTS, CHB, S0, S1 ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), & VCOSMO2SOLVENT(NUMPOINTS), SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS),& TM(NUMPOINTS), DELTAHF(NUMPOINTS), DELTACP(NUMPOINTS), SYSTEMP(NUMPOINTS), & X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS), SOLUTEFILE(NUMPOINTS), SOLVENT1FILE(NUMPOINTS), & UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS), SIGMASOLUTE_NHB(COMPSEG,NUMPOINTS), & SIGMASOLUTE_HB(COMPSEG,NUMPOINTS), COUNTERSOLUTE(COMPSEG,NUMPOINTS), ACOSMOSOLUTE(NUMPOINTS), & SIGMA1SOLVENT_NHB(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_HB(COMPSEG,NUMPOINTS), & COUNTER1SOLVENT(COMPSEG,NUMPOINTS), SIGMA2SOLVENT_HB(COMPSEG,NUMPOINTS), & SIGMA2SOLVENT_NHB(COMPSEG,NUMPOINTS), COUNTER2SOLVENT(COMPSEG,NUMPOINTS), & ACOSMO2SOLVENT(NUMPOINTS), SOLV2(NUMPOINTS), & UNIT2SOLV(NUMPOINTS), X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), & SOLVENT2FILE(NUMPOINTS), STAT = ISTAT) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMTEMP, NUMPOINTS, CHB, S0, S1 DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), SOLV2(L), & VCOSMO2SOLVENT(L), TM(L), DELTAHF(L), DELTACP(L), SYSTEMP(L), X2DATA(L)!, X3DATA(L) END DO CLOSE (12) DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNIT2SOLV (I) = I + 1020 UNITSOLU (I) = I + 10020 SOLVENT1FILE (I) = TRIM(SOLV1(I))//'_refined.txt' ![USER SPECIFIED FILE LOCATION] SOLVENT2FILE (I) = TRIM(SOLV2(I))//'_refined.txt' ![USER SPECIFIED FILE LOCATION] SOLUTEFILE (I) = TRIM(SOLU(I))//'_refined.txt' ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLVENT2FILE(I), SOLUTEFILE(I), X2DATA(I), TM(I) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT_NHB(J,K), SIGMA1SOLVENT_HB(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT_NHB(J,K) + SIGMA1SOLVENT_HB(J,K) END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT2SOLV(K), FILE = SOLVENT2FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO2SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT2SOLV(K),*) COUNTER2SOLVENT(J,K), SIGMA2SOLVENT_NHB(J,K), SIGMA2SOLVENT_HB(J,K) ACOSMO2SOLVENT(K) = ACOSMO2SOLVENT(K) + SIGMA2SOLVENT_NHB(J,K) + SIGMA2SOLVENT_HB(J,K) END DO CLOSE (UNIT2SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE_NHB(J,K), SIGMASOLUTE_HB(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE_NHB(J,K) + SIGMASOLUTE_HB(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins_refined.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_refined_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_refined_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A7,1X,A7) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)", "COMP(3)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","X3" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP, 2), COUNTER(COMPSEG),& MIXPROFILE(COMPSEG,2), DELTAW(COMPSEG,COMPSEG,2,2),& SEGGAMMA(COMPSEG,2),SEGGAMMAOLD(COMPSEG,2),CONVERG(COMPSEG,2),SEGGAMMAPURE(COMPSEG,COMP,2),& SEGGAMMAPUREOLD(COMPSEG,COMP,2), CONVPURE(COMPSEG,COMP), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP, 2), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS ! BEGIN ITERATION LOOP FOR ALL POINTS (température), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMTEMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 DO I = 1, COMPSEG COUNTER(I) = COUNTER1SOLVENT(I,P) !COUNTER(I) = COUNTER1SOLVENT(I,1) ! 1->NHB, 2->HB !COUNTER(I) = COUNTER1SOLVENT(I,2) ! SIGMA(segment, species, nhb or hb) SIGMA(I,1,1) = SIGMA1SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S1**2)))& + SIGMA1SOLVENT_NHB(I,P) SIGMA(I,2,1) = SIGMA2SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S1**2)))& + SIGMA2SOLVENT_NHB(I,P) SIGMA(I,3,1) = SIGMASOLUTE_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMASOLUTE_NHB(I,P) SIGMA(I,1,2) = SIGMA1SOLVENT_HB(I,P)*(1-dexp(-COUNTER(I)**2/(2*S1**2))) SIGMA(I,2,2) = SIGMA2SOLVENT_HB(I,P)*(1-dexp(-COUNTER(I)**2/(2*S1**2))) SIGMA(I,3,2) = SIGMASOLUTE_HB(I,P)*(1-dexp(-COUNTER(I)**2/(2*S0**2))) !write(*,*) (dexp(-COUNTER(I)**2/(2*S0**2))) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMO2SOLVENT(P) ACOSMO(3) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMO2SOLVENT(P) VCOSMO(3) = VCOSMOSOLUTE(P) !X(1) = 1.0D0 - X2DATA(P) - X3DATA(P) !X(2) = X2DATA(P) X(2) = X2DATA(P) !1.0D0 - X2DATA(P) - X3DATA(P) X(1) = 1 - X(2) X(3) = 0.000000001D0 !X3DATA(P) !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) - X(3) X2OLD = X(2) X3OLD = X(3) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] ! MIXPROFILE IS p(sigmai(J)) FOR THE MIXTURE in case NHB (1), HB (2) MIXPROFILE(J,1) = ((X(1)*SIGMA(J,1,1))+(X(2)*SIGMA(J,2,1))+(X(3)*SIGMA(J,3,1))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) MIXPROFILE(J,2) = ((X(1)*SIGMA(J,1,2))+(X(2)*SIGMA(J,2,2))+(X(3)*SIGMA(J,3,2))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO H=1,2 DO B=1,2 DO I = 1, COMPSEG DO K = 1, COMPSEG IF (H==B.AND.H==2.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN ! DELTAW(I,K,H,B) is delta for segment I for H=1 (NHB) or H=2 (HB) and segment K DELTAW(I,K,H,B) = (ALPHAPRIME/2.0D0)*(COUNTER(I)+COUNTER(K))**2.0D0 - & CHB*(COUNTER(I)-COUNTER(K))**2 ELSE DELTAW(I,K,H,B) = (ALPHAPRIME/2.0D0)*(COUNTER(I)+COUNTER(K))**2.0D0 ENDIF END DO END DO ENDDO ENDDO ! ITERATION FOR MIXTURE SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 SEGGAMMA = 1.0D0 !INTIAL VALUE FOR MIXTURE SEGMENT ACTIVITY COEF. DO SEGGAMMAOLD = SEGGAMMA DO B=1,2 DO I = 1, COMPSEG SUMMATION = 0.0D0 DO H=1,2 DO K = 1, COMPSEG SUMMATION = SUMMATION + MIXPROFILE(K,H)*SEGGAMMAOLD(K,H) * & DEXP(-DELTAW(I,K,B,H)/(RGAS * SYSTEMP(P))) ENDDO ENDDO SEGGAMMA(I,B) = DEXP(-DLOG(SUMMATION)) SEGGAMMA(I,B) = (SEGGAMMA(I,B) + SEGGAMMAOLD(I,B))/2.0D0 END DO ENDDO DO B=1, 2 DO I=1, COMPSEG CONVERG(I,B) = DABS((SEGGAMMA(I,B) - SEGGAMMAOLD(I,B))/SEGGAMMAOLD(I,B)) END DO ENDDO IF (MAXVAL(CONVERG) <=0.000001) EXIT ENDDO ! ITERATION FOR PURE SPECIES, L, SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 ! ITERATION FOR SEGMENT ACITIVITY COEF (PURE SPECIES) ! HERE p(sigma) is just SIGMA/Ai total IF (ITER==1) THEN DO L = 1, COMP SEGGAMMAPURE (:,L,:) = 1.0D0 DO SEGGAMMAPUREOLD (:,L,:) = SEGGAMMAPURE (:,L,:) DO I = 1, COMPSEG DO B=1,2 SUMMATION = 0.0D0 DO H=1,2 DO K = 1, COMPSEG SUMMATION = SUMMATION + (SIGMA(K,L,H)/ACOSMO(L)) * & SEGGAMMAPUREOLD(K,L,H)*DEXP(-DELTAW(I,K,B,H)/(RGAS*SYSTEMP(P))) ENDDO END DO SEGGAMMAPURE(I,L,B)=DEXP(-DLOG(SUMMATION)) SEGGAMMAPURE(I,L,B)=(SEGGAMMAPURE(I,L,B)+SEGGAMMAPUREOLD(I,L,B))/2.0D0 ENDDO END DO DO B=1,2 DO I=1,COMPSEG CONPR(I,L,B)=DABS((SEGGAMMAPURE(I,L,B)-SEGGAMMAPUREOLD(I,L,B))/& SEGGAMMAPUREOLD(I,L,B)) ENDDO ENDDO IF (MAXVAL(CONPR) <=0.000001) EXIT ENDDO ENDDO END IF ! THE STAVERMAN-GUGGENHEIM EQUATION DO I = 1,COMP RNORM(I) = VCOSMO(I)/VNORM QNORM(I) = ACOSMO(I)/ANORM END DO DO I = 1, COMP THETA(I) = (X(I)*QNORM(I))/(X(1)*QNORM(1) + X(2)*QNORM(2) + X(3)*QNORM(3)) PHI(I) = (X(I)*RNORM(I))/(X(1)*RNORM(1) + X(2)*RNORM(2) + X(3)*RNORM(3)) LSG(I) = (COORD/2.0D0)*(RNORM(I)-QNORM(I))-(RNORM(I)-1.0D0) END DO ! GAMMASG1 AND GAMMASG2 ARE ACTUALLY LNGAMMASG DO I = 1, COMP LNGAMMASG(I) = DLOG(PHI(I)/X(I)) + (COORD/2.0D0)*QNORM(I)* & DLOG(THETA(I)/PHI(I)) + LSG(I) - (PHI(I)/X(I))*(X(1)*LSG(1) + X(2)*LSG(2) + X(3)*LSG(3)) END DO !CALCULATION OF GAMMAS SUMGAMMA = 0.0D0 DO K = 1, COMP DO B = 1, 2 DO I=1, COMPSEG SUMGAMMA(K) = SUMGAMMA(K)+((SIGMA(I,K,B)/AEFFPRIME)*(DLOG(SEGGAMMA(I,B)/ & SEGGAMMAPURE(I,K,B)))) END DO ENDDO GAMMA(K) = DEXP(SUMGAMMA(K) + (LNGAMMASG(K))) LNGAMMA(K) = DLOG(GAMMA(K)) END DO ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(3) = DEXP((DELTAHF(P)/8.3145D0/TM(P))*(1.0D0-TM(P)/SYSTEMP(P)) - & (DELTACP(P)/8.3145D0)*(DLOG(TM(P)/SYSTEMP(P))-(TM(P)/SYSTEMP(P))+1.0D0) & - LNGAMMA(3)) !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(3) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(3) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(3) = (X(3) + 2.0D0*X3OLD)/3.0D0 MFSUM = X(1)+X(2)+X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM X(3) = X(3)/MFSUM ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(3) - X3OLD)/(X3OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),X(3) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(3) write (16,*) GAMMA(3) write(*,*) X(3) ! #### CAS OU ON INCREMENTE LA TEMPERATURE : SYSTEMP(P) = SYSTEMP(P) + 5 END DO ! Iteration temperature (variable P) END DO ! Iteration points (variable R) CLOSE (13) ! CLOSE OUTPUT FILE CLOSE (15) ! CLOSE l'AUTRE OUTPUT FILE CLOSE (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, DELTACP, SYSTEMP, & X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE_NHB, SIGMASOLUTE_HB, COUNTERSOLUTE, ACOSMOSOLUTE,& SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_HB,COUNTER1SOLVENT, & SIGMA2SOLVENT_NHB, SIGMA2SOLVENT_HB, VCOSMO2SOLVENT, & SOLV2, SOLVENT2FILE, ACOSMO2SOLVENT, UNIT2SOLV, X3DATA, & ACOSMO1SOLVENT, X1DATA, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" !SOLV2, VCOSMO2SOLVENT, SOLVENT2FILE, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT !X1DATA, X3DATA, UNIT2SOLV, END PROGRAM GAMMASOLUBILITY I.3/COSMO-SAC (2010) I.3.1/COSMO-SAC 2010 binaire I.3.1.1/COSMO-SAC 2010 binaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PHYSICAL CONSTANTS AND PARAMETERS: ! EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom) ! S0 = SHAPE PARAMETER ? ! 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) ! 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 ! VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3) ! ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM ! OF THE INDIVIDUAL PROFILE. ! ! IN DATA.txt FILE !!!!!!!!! ! ! PROGRAM CURRENTLY SETUP FOR PURE SOLVENT SOLUBILITY PREDICTIONS ONLY !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.25D0 ! ou 7.5 selon publi ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 !REAL*8, PARAMETER :: S0 = 0.05D0 !0.007D0 !0.007D0 ! shape parameter REAL*8 :: S0, S1 REAL*8, PARAMETER :: AES = 6525.69D0 ! MODIFICATION OF LIN CES REAL*8, PARAMETER :: BES = 148590000.0D0 ! MODIFICATION OF LIN CES !REAL*8, PARAMETER :: COH_OH = 0.0D0! 4013.78D0 ! OH-OH INTERACTION PARAMETER !REAL*8, PARAMETER :: COT_OT = 0.0D0! 932.31D0 ! OT-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OT = 0.0D0! 3016.43D0 ! OH-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OH = 4013.78D0 ! OH-OH INTERACTION PARAMETER !REAL*8, PARAMETER :: COT_OT = 932.31D0 ! OT-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OT = 3016.43D0 ! OH-OT INTERACTION PARAMETER REAL*8 :: COH_OH REAL*8 :: COT_OT REAL*8 :: COH_OT INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:,:,:), ALLOCATABLE :: DELTAW REAL*8, DIMENSION (:,:,:), ALLOCATABLE :: SIGMA, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONPR REAL*8, DIMENSION (:,:), ALLOCATABLE :: WHB, SEGGAMMA, SEGGAMMAOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE_NHB, SIGMASOLUTE_OH, SIGMASOLUTE_OT, SIGMASOLUTE_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_OH, SIGMA1SOLVENT_OT, SIGMA1SOLVENT_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: MIXPROFILE, COUNTERSOLUTE REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONVERG REAL*8, DIMENSION (:), ALLOCATABLE :: X1DATA, X2DATA, X, COUNTER !, X3DATA, REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, DELTAHF, DELTACP REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, SUMMATION , CES REAL*8 :: RELTOL, CHANGE, X2OLD, X1NEW, X2NEW, MFSUM ! X3NEW INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMTEMP, H, B INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU ! UNIT2SOLV, CHARACTER (40) :: INPUTFILE, OUTPUTFILE CHARACTER (30), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE !, SOLVENT2FILE, CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLV1 !, SOLV2 CHARACTER (20), DIMENSION (:), ALLOCATABLE :: SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 !FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS FPOL = 0.6916D0 !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 2 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data_10refined_cp.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMTEMP = 41 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMTEMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 CLOSE(12) write(*,*) NUMTEMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), VCOSMO2SOLVENT(NUMPOINTS),& SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS), TM(NUMPOINTS), DELTAHF(NUMPOINTS),& SYSTEMP(NUMPOINTS), X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS),SOLUTEFILE(NUMPOINTS),& SOLVENT1FILE(NUMPOINTS), UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS),& SIGMASOLUTE_NHB(COMPSEG,NUMPOINTS), SIGMASOLUTE_OH(COMPSEG,NUMPOINTS),& SIGMASOLUTE_OT(COMPSEG,NUMPOINTS), SIGMASOLUTE_HB(COMPSEG,NUMPOINTS),& COUNTERSOLUTE(COMPSEG,NUMPOINTS), ACOSMOSOLUTE(NUMPOINTS), DELTACP(NUMPOINTS),& SIGMA1SOLVENT_NHB(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_OH(COMPSEG,NUMPOINTS), & SIGMA1SOLVENT_OT(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_HB(COMPSEG,NUMPOINTS), & COUNTER1SOLVENT(COMPSEG,NUMPOINTS), STAT = ISTAT) !, & !SIGMA2SOLVENT(COMPSEG,NUMPOINTS), COUNTER2SOLVENT(COMPSEG,NUMPOINTS), & !ACOSMO2SOLVENT(NUMPOINTS)), SOLV2(NUMPOINTS), !UNIT2SOLV(NUMPOINTS),X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), SOLVENT2FILE(NUMPOINTS)) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMTEMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), TM(L), DELTAHF(L), & SYSTEMP(L), DELTACP(L), X2DATA(L) ! , X1DATA(L), X3DATA(L), SOLV2(L), VCOSMO2SOLVENT(L) END DO CLOSE (12) DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNITSOLU (I) = I + 10020 SOLVENT1FILE (I) = TRIM(SOLV1(I))//'_10refined.txt' ![USER SPECIFIED FILE LOCATION] SOLUTEFILE (I) = TRIM(SOLU(I))//'_10refined.txt' ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLUTEFILE(I), X2DATA(I), TM(I), DELTAHF(I), SYSTEMP(I), VCOSMOSOLUTE(I) END DO ! SIGMA1SOLVENT_HB = 0.0D0 ! SIGMASOLUTE_HB = 0.0D0 DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT_NHB(J,K), & SIGMA1SOLVENT_OH(J,K), SIGMA1SOLVENT_OT(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT_NHB(J,K) + & SIGMA1SOLVENT_OH(J,K) + SIGMA1SOLVENT_OT(J,K) SIGMA1SOLVENT_HB(J,K) = SIGMA1SOLVENT_OH(J,K) + SIGMA1SOLVENT_OT(J,K) END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE_NHB(J,K),& SIGMASOLUTE_OH(J,K), SIGMASOLUTE_OT(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE_NHB(J,K) +& SIGMASOLUTE_OH(J,K) + SIGMASOLUTE_OT(J,K) SIGMASOLUTE_HB(J,K) = SIGMASOLUTE_OH(J,K) + SIGMASOLUTE_OT(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins_10refined.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_10refined_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_10refined_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "control.txt" OPEN (UNIT = 17, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A7,1X,A7) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","GAMMA1","GAMMA2" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP, 3), COUNTER(COMPSEG),& MIXPROFILE(COMPSEG,3), DELTAW(COMPSEG,COMPSEG,3,3),& SEGGAMMA(COMPSEG,3),SEGGAMMAOLD(COMPSEG,3),CONVERG(COMPSEG,3),SEGGAMMAPURE(COMPSEG,COMP,3),& SEGGAMMAPUREOLD(COMPSEG,COMP,3), CONVPURE(COMPSEG,COMP), WHB(COMPSEG, COMPSEG), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP, 3), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS ! write(*,*) SYSTEMP(P) ! BEGIN ITERATION LOOP FOR ALL POINTS (température), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMTEMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 WHB = 0.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 ! EXPRESSION OF CES CES = AES + BES/(SYSTEMP(P)**2) !CES = ALPHAPRIME/2 DO I = 1, COMPSEG !COUNTER(I) = COUNTER1SOLVENT(I,P) COUNTER(I) = COUNTER1SOLVENT(I,P) ! 1->NHB, 2->OH, 3->OT ! SIGMA(SEGMENT, SPECIES, NHB or OH or OT) SIGMA(I,1,1) = SIGMA1SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMA1SOLVENT_NHB(I,P) SIGMA(I,2,1) = SIGMASOLUTE_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S1**2)))& + SIGMASOLUTE_NHB(I,P) SIGMA(I,1,2) = SIGMA1SOLVENT_OH(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,2,2) = SIGMASOLUTE_OH(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S1**2))) SIGMA(I,1,3) = SIGMA1SOLVENT_OT(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,2,3) = SIGMASOLUTE_OT(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S1**2))) ! SIGMA(I,1,1) = SIGMA1SOLVENT_NHB(I,P) ! SIGMA(I,2,1) = SIGMASOLUTE_NHB(I,P) ! SIGMA(I,1,2) = SIGMA1SOLVENT_OH(I,P) ! SIGMA(I,2,2) = SIGMASOLUTE_OH(I,P) ! SIGMA(I,1,3) = SIGMA1SOLVENT_OT(I,P) ! SIGMA(I,2,3) = SIGMASOLUTE_OT(I,P) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMOSOLUTE(P) X(1) = 1.0D0 - X2DATA(P) X(2) = X2DATA(P) !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) X2OLD = X(2) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] ! MIXPROFILE IS p(sigmai(J)) FOR THE MIXTURE in case NHB (1), OH (2), OT(3) MIXPROFILE(J,1) = ((X(1)*SIGMA(J,1,1))+(X(2)*SIGMA(J,2,1))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) MIXPROFILE(J,2) = ((X(1)*SIGMA(J,1,2))+(X(2)*SIGMA(J,2,2))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) MIXPROFILE(J,3) = ((X(1)*SIGMA(J,1,3))+(X(2)*SIGMA(J,2,3))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO H=1,3 DO B=1,3 DO I = 1, COMPSEG DO K = 1, COMPSEG IF (H==B.AND.H==2.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN ! DELTAW(I,K,H,B) is delta for segment I for H=1 (NHB) or H=2 (OH) H=3 (OT) for segment K DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COH_OH*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (H==B.AND.H==3.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COT_OT*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (((H==2.AND.B==3).OR.(H==3.AND.B==2)).AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COH_OT*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (H==1.OR.B==1.OR.(COUNTER(I)*COUNTER(K))>=0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 ENDIF END DO END DO ENDDO ENDDO ! ITERATION FOR MIXTURE SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 SEGGAMMA = 1.0D0 !INTIAL VALUE FOR MIXTURE SEGMENT ACTIVITY COEF. DO SEGGAMMAOLD = SEGGAMMA DO B=1,3 DO I = 1, COMPSEG SUMMATION = 0.0D0 DO H=1,3 DO K = 1, COMPSEG SUMMATION = SUMMATION + MIXPROFILE(K,H)*SEGGAMMAOLD(K,H) * & DEXP(-DELTAW(I,K,B,H)/(RGAS * SYSTEMP(P))) ENDDO ENDDO SEGGAMMA(I,B) = DEXP(-DLOG(SUMMATION)) SEGGAMMA(I,B) = (SEGGAMMA(I,B) + SEGGAMMAOLD(I,B))/2.0D0 END DO ENDDO DO B=1, 3 DO I=1, COMPSEG CONVERG(I,B) = DABS((SEGGAMMA(I,B) - SEGGAMMAOLD(I,B))/SEGGAMMAOLD(I,B)) END DO ENDDO IF (MAXVAL(CONVERG) <=0.000001) EXIT ENDDO ! ITERATION FOR PURE SPECIES, L, SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 ! ITERATION FOR SEGMENT ACITIVITY COEF (PURE SPECIES) ! HERE p(sigma) is just SIGMA/Ai total IF (ITER==1) THEN DO L = 1, COMP !SEGGAMMAPURE (:,L,:) = 1.0D0 SEGGAMMAPURE = 1.0D0 DO !SEGGAMMAPUREOLD (:,L,:) = SEGGAMMAPURE (:,L,:) SEGGAMMAPUREOLD = SEGGAMMAPURE DO I = 1, COMPSEG DO B=1,3 SUMMATION = 0.0D0 DO H=1,3 DO K = 1, COMPSEG SUMMATION = SUMMATION + (SIGMA(K,L,H)/ACOSMO(L)) * & SEGGAMMAPUREOLD(K,L,H)*DEXP(-DELTAW(I,K,B,H)/(RGAS*SYSTEMP(P))) ENDDO END DO SEGGAMMAPURE(I,L,B)=DEXP(-DLOG(SUMMATION)) SEGGAMMAPURE(I,L,B)=(SEGGAMMAPURE(I,L,B)+SEGGAMMAPUREOLD(I,L,B))/2.0D0 ENDDO END DO DO B=1,3 DO I=1,COMPSEG CONPR(I,L,B)=DABS((SEGGAMMAPURE(I,L,B)-SEGGAMMAPUREOLD(I,L,B))/& SEGGAMMAPUREOLD(I,L,B)) ENDDO ENDDO IF (MAXVAL(CONPR) <=0.000001) EXIT ENDDO ENDDO ENDIF ! THE STAVERMAN-GUGGENHEIM EQUATION DO I = 1,COMP RNORM(I) = VCOSMO(I)/VNORM QNORM(I) = ACOSMO(I)/ANORM END DO DO I = 1, COMP THETA(I) = (X(I)*QNORM(I))/(X(1)*QNORM(1) + X(2)*QNORM(2)) PHI(I) = (X(I)*RNORM(I))/(X(1)*RNORM(1) + X(2)*RNORM(2)) LSG(I) = (COORD/2.0D0)*(RNORM(I)-QNORM(I))-(RNORM(I)-1.0D0) END DO ! GAMMASG1 AND GAMMASG2 ARE ACTUALLY LNGAMMASG DO I = 1, COMP LNGAMMASG(I) = DLOG(PHI(I)/X(I)) + (COORD/2.0D0)*QNORM(I)* & DLOG(THETA(I)/PHI(I)) + LSG(I) - (PHI(I)/X(I))*(X(1)*LSG(1) + X(2)*LSG(2)) END DO !CALCULATION OF GAMMAS SUMGAMMA = 0.0D0 DO K = 1, COMP DO B = 1, 3 DO I=1, COMPSEG SUMGAMMA(K) = SUMGAMMA(K)+((SIGMA(I,K,B)/AEFFPRIME)*(DLOG(SEGGAMMA(I,B)/ & SEGGAMMAPURE(I,K,B)))) END DO ENDDO GAMMA(K) = DEXP(SUMGAMMA(K) + (LNGAMMASG(K))) LNGAMMA(K) = DLOG(GAMMA(K)) END DO ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(2) = DEXP((DELTAHF(P)/8.3145D0/TM(P))*(1.0D0-TM(P)/SYSTEMP(P))- & (DELTACP(P)/8.3145D0)*(DLOG(TM(P)/SYSTEMP(P))-(TM(P)/SYSTEMP(P))+1.0D0) & - LNGAMMA(2)) ! Affichage des résultats intermédiaires ! WRITE(*,*) P, X(2) !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(2) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(2) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(2) = (X(2) + 2.0D0*X2OLD)/3.0D0 MFSUM = X(1)+X(2) ! +X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM ! Pour afficher la fraction molaire finale intermédiaire ! write(*,*) 'final', X(2) ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(2) - X2OLD)/(X2OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 write(17,*) X2OLD, X(2), LNGAMMASG(2), SUMGAMMA(1) !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),GAMMA(1),GAMMA(2) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(2) write (16,*) GAMMA(2) write(*,*) X(2) ! #### CAS OU ON INCREMENTE LA TEMPERATURE : SYSTEMP(P) = SYSTEMP(P) + 5 END DO ! Iteration temperature (variable P) END DO ! Iteration points (variable R) CLOSE (13) ! CLOSE OUTPUT FILE close (15) ! CLOSE l'AUTRE OUTPUT FILE close (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, WHB, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, SYSTEMP, DELTACP,& X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE_NHB, SIGMASOLUTE_OH, SIGMASOLUTE_OT, SIGMASOLUTE_HB, & COUNTERSOLUTE, ACOSMOSOLUTE, SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_OH, & COUNTER1SOLVENT, SIGMA1SOLVENT_OT, SIGMA1SOLVENT_HB, ACOSMO1SOLVENT, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" !SOLV2, VCOSMO2SOLVENT, SOLVENT2FILE, SIGMA2SOLVENT, COUNTER2SOLVENT, ACOSMO2SOLVENT !X1DATA, X3DATA, UNIT2SOLV, END PROGRAM GAMMASOLUBILITY I.3.1.2/COSMO-SAC 2010 binaire MATLAB function [y] = cosmosac10(X) %[coh-oh cotot cohot S0 S1]) % % Cosmo-sac renvoie l'erreur entre la solubilité attendue et celle trouvée % on rentre comme parametre : Cohho Cotot Cohot S0 et S1 dans le but de les optimiser %format long; % coh=X(1); %interaction OH_OH cot=X(2); %interaction OT_OT coh_ot=X(3); %interaction OH_OT S0=X(4); %probabilité HB pour solvant S1=X(5); %probabilité HB pour soluté AES= 8806.9; %original 6525.69; BES= 7.5852E7; %original 1.4859E8; comp=2; compseg = 51;% NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) EO = 0.0002395 ;% PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) R = 0.001987 ;% KCAL/(MOL K) aeffprime = 7.2500 ;% ANGSTROMS^2 ANORM = 79.5300 ;% NORMALIZED SURFACE AREA, ANGSTROMS^2 VNORM = 66.6900 ;% NORMALIZED CAVITY VOLUME, ANGSTROMS^3 COORD = 10.0000 ;% COORDINATE NUMBER, Z, KLAMT SET TO 7.2 % % Lecture des données d'entrée et initialisation de la boucle A = importdata('input2.txt'); % A(p,1)=solu,2=solv,3=Vcosmosolu,4=vcosmosolv,5=Tm,6=Hm,7=T,8=X2 % A.textdata{p,1}=solu.txt, (2)=solv.txt %A.data(p,1)=Vcosmosolu, (p,2)=Vcosmosolv, (p,3)=Tm, (p,4)=Hm, (p,5)=T, (p,6)=X2 N = size(A.data,1) sol=zeros(N,1); solub=zeros(N,1); mixprofile=zeros(compseg,2); x=zeros(2,1); sigmasol=zeros(compseg,2,3); retol = 1E-6; %sigmasolv=zeros(compseg,2,2); % Solute = load(A.textdata{1,1}); Solv = load(A.textdata{1,2}); for p=1:N CES=AES+BES/(A.data(p,5)^2); % ajouter une condition sur le cas précédent pour pas tout refaire... % w=p-1; % if (p~=1&&A.textdata{p,1}~=A.textdata{w,1}) %aime pas le (p-1) en argument... Solute = load(A.textdata{p,1}); % end Solv = load(A.textdata{p,2}); Delta=zeros(51,51,2,2); Acosmo=zeros(2); Acosmo(1)=ones(1,51)*Solv(:,2)+ones(1,51)*Solv(:,3)+ones(1,51)*Solv(:,4); Acosmo(2)=ones(1,51)*Solute(:,2)+ones(1,51)*Solute(:,3)+ones(1,51)*Solute(:,4); sigmasol(:,1,1)=Solv(:,2)+(Solv(:,3)+Solv(:,4)).*(exp(-(Solv(:,1).^2)./(2*S0^2))); sigmasol(:,2,1)=Solute(:,2)+(Solute(:,3)+Solute(:,4)).*(exp(-(Solute(:,1).^2)./(2*S1^2))); sigmasol(:,1,2)=Solv(:,3).*(1-exp(-(Solv(:,1).^2)./(2*S0^2))); sigmasol(:,2,2)=Solute(:,3).*(1-exp(-(Solute(:,1).^2)./(2*S1^2))); sigmasol(:,1,3)=Solv(:,4).*(1-exp(-(Solv(:,1).^2)./(2*S0^2))); sigmasol(:,2,3)=Solute(:,4).*(1-exp(-(Solute(:,1).^2)./(2*S1^2))); x(2)=A.data(p,6); x(1)=1-x(2); iter=-1; %while err > 10^-8 % boucle sur l'erreur entre solubilité prédite et exp change = 1; while change > retol % boucle sur la détermination de x(2) iter = iter +1 ; % Calcul du gamma combinatoire lngammasg = lngammac(A.data(p,1),A.data(p,2),Acosmo(1),Acosmo(2),x(2),VNORM,ANORM,COORD); % % Calcul du gamma résiduel mixprofile(:,1)=((x(1)*sigmasol(:,1,1))+(x(2)*sigmasol(:,2,1)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); mixprofile(:,2)=((x(1)*sigmasol(:,1,2))+(x(2)*sigmasol(:,2,2)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); mixprofile(:,3)=((x(1)*sigmasol(:,1,3))+(x(2)*sigmasol(:,2,3)))./((x(1)*Acosmo(1)+x(2)*Acosmo(2))); %x(1)=1-x(2); x2old=x(2); for h=1:3 for b=1:3 for i=1:compseg for k=1:compseg Delta(i,k,h,b)=CES*(Solv(i,1)+Solv(k,1))^2; if (h==b&&h==2&&(Solv(i,1)*Solv(k,1))<0) Delta(i,k,h,b)=CES*(Solv(i,1)+Solv(k,1))^2-coh*(Solv(i,1)-Solv(k,1))^2; end if (h==b&&h==3&&(Solv(i,1)*Solv(k,1))<0) Delta(i,k,h,b)=CES*(Solv(i,1)+Solv(k,1))^2-cot*(Solv(i,1)-Solv(k,1))^2; end if (((h==3&&b==2)|(h==2&&b==3))&&(Solv(i,1)*Solv(k,1))<0) Delta(i,k,h,b)=CES*(Solv(i,1)+Solv(k,1))^2-coh_ot*(Solv(i,1)-Solv(k,1))^2; end end end end end seggamma = ones(compseg,3); converg=ones(compseg,3); while max(max(converg))>=0.00001 % boucle pour seggamma du mélange seggammaold = seggamma; for b=1:3 for i=1:compseg summation = 0; for h=1:3 for k=1:compseg summation = summation + mixprofile(k,h)*seggammaold(k,h)*... exp(-Delta(i,k,b,h)/(R*A.data(p,5))); end % en retirant la boucle sur k: summation = summation + mixprofile(:,h)'*(seggammaold(:,h).*exp(-Delta(i,:,b,h)'./(R*A.data(p,5)))); % malheureusement, c'est moins intéressant niveau optimisation. pkoi ??? end seggamma(i,b)=exp(-log(summation)); seggamma(i,b)=(seggamma(i,b)+2*seggammaold(i,b))/3; end end converg=abs((seggamma-seggammaold)./seggammaold); end % Fin boucle % if iter==0 conpr=ones(51,comp,3); seggammapureold=ones(compseg,comp,3); for l=1:comp seggammapure=ones(compseg,comp,3); while (max(max(max(conpr(:,l,:))))>=0.00001) % début boucle pour calcul de seggammapure conpr(:,l,:)=zeros(51,3); %zeros(51,2,2); seggammapureold(:,l,:)=seggammapure(:,l,:); for i=1:compseg for b=1:3 summation=0; for h=1:3 for k=1:51 summation=summation+(sigmasol(k,l,h)/Acosmo(l))... *seggammapureold(k,l,h)*exp(-Delta(i,k,b,h)/(R*A.data(p,5))); end end seggammapure(i,l,b)=exp(-log(summation)); seggammapure(i,l,b)=(seggammapure(i,l,b)+seggammapureold(i,l,b))/2; end end conpr(:,l,:)=abs((seggammapure(:,l,:)-seggammapureold(:,l,:))./seggammapureold(:,l,:)); end %fin boucle seggammapure (WHILE) end %fin boucle sur l end % % Calcul de gamma residuel gamma=ones(2,1); lngamma=ones(2,1); sumgamma=zeros(2,1); for k=1:comp for b=1:3 for i=1:compseg sumgamma(k)=sumgamma(k) + ((sigmasol(i,k,b)/aeffprime)*(log(seggamma(i,b)/... seggammapure(i,k,b)))); end end gamma(k)=exp(sumgamma(k)+(lngammasg(k))); lngamma(k)=log(gamma(k)); end % Recalcul de la fraction molaire x(2) = exp((A.data(p,4)/8.3145/A.data(p,3))*(1-A.data(p,3)/A.data(p,5))-lngamma(2)); x(2) = (x(2)+2*x2old)/3; mfsum=x(1)+x(2); x(1)=x(1)/mfsum; x(2)=x(2)/mfsum; change=abs((x(2)-x2old)/(x2old+1E-16)); end sol(p)=x(2); solub(p)=(log(A.data(p,6))-log(sol(p)))/sqrt(N); end %y=sol; sol; disp(X); y=solub; end function [y]=lngammac(Vcosmosolute,Vcosmosolv,Acosmosolv,Acosmosolute,Xs,Vnorm,Anorm,COORD) %theta=zeros(2,1); %phi=zeros(2,1); %lsg=zeros(2,1); %lngammasg=zeros(2,1); X(1)=1-Xs; X(2)=Xs; Rnorm(1)=Vcosmosolv/Vnorm; Rnorm(2)=Vcosmosolute/Vnorm; Qnorm(1)=Acosmosolv/Anorm; Qnorm(2)=Acosmosolute/Anorm; theta=(X.*Qnorm)./(X*Qnorm'); phi=(X.*Rnorm)./(X*Rnorm'); lsg=(COORD/2).*(Rnorm-Qnorm)-(Rnorm-1); lngammasg=log(phi./X)+(COORD/2).*Qnorm.*log(theta./phi)+lsg... - (phi./X)*(X*lsg'); %[X(2),Rnorm(2), Qnorm(2), theta(2), phi(2), lsg(2)] y=lngammasg; end I.3.2/COSMO-SAC 2010 ternaire FORTRAN PROGRAM GAMMASOLUBILITY !*********************************************************************************** ! This program, GAMMASOLUBILITY, calculates activity coefficients for a solute and ! solvent(s) for pure solvent systems and binary or mixed solvent systems. The ! activity coefficients are then used to predict a solubility value for the solute ! in that mixture. ! ! This program uses the COSMO-SAC model as published (Lin, S.T., ! S.I. Sandler, Ind. Eng. Chem. Res. 41, (2002), 899-913). ! ! PHYSICAL CONSTANTS AND PARAMETERS: ! EO = PERMITTIVITY IN A VACUUM (e**2*mol/Kcal*Angstrom) ! S0 = SHAPE PARAMETER ? ! 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) ! 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 ! VCOSMO = CAVITY VOLUME FROM COSMO OUTPUT (A**3) ! ACOSMO = MOLECULAR SURFACE AREA FROM COSMO OUTPUT (A**2) --THE SUM ! OF THE INDIVIDUAL PROFILE. ! ! IN DATA.txt FILE !!!!!!!!! ! ! PROGRAM CURRENTLY SETUP FOR PURE SOLVENT SOLUBILITY PREDICTIONS ONLY !********************************************************************************** IMPLICIT NONE REAL*8, PARAMETER :: EO = 0.0002395D0 ! PERMITTIVITY OF FREE SPACE (e^2*MOL/KCAL/ANG^2) REAL*8, PARAMETER :: AEFFPRIME = 7.25D0 ! ou 7.5 selon publi ANGSTROMS^2 REAL*8, PARAMETER :: RGAS = 0.001987D0 ! KCAL/(MOL K) REAL*8, PARAMETER :: VNORM = 66.69D0 ! NORMALIZED CAVITY VOLUME, ANGSTROMS^3 REAL*8, PARAMETER :: ANORM = 79.53D0 ! NORMALIZED SURFACE AREA, ANGSTROMS^2 REAL*8, PARAMETER :: COORD = 10.0D0 ! COORDINATE NUMBER, Z, KLAMT SET TO 7.2 REAL*8, PARAMETER :: AES = 6525.69D0 ! MODIFICATION OF LIN CES REAL*8, PARAMETER :: BES = 148590000.0D0 ! MODIFICATION OF LIN CES !REAL*8, PARAMETER :: COH_OH = 0.0D0! 4013.78D0 ! OH-OH INTERACTION PARAMETER !REAL*8, PARAMETER :: COT_OT = 0.0D0! 932.31D0 ! OT-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OT = 0.0D0! 3016.43D0 ! OH-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OH = 4013.78D0 ! OH-OH INTERACTION PARAMETER !REAL*8, PARAMETER :: COT_OT = 932.31D0 ! OT-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: COH_OT = 3016.43D0 ! OH-OT INTERACTION PARAMETER REAL*8 :: COH_OH ! OH-OH INTERACTION PARAMETER REAL*8 :: COT_OT ! OT-OT INTERACTION PARAMETER REAL*8 :: COH_OT ! OH-OT INTERACTION PARAMETER !REAL*8, PARAMETER :: S0 = 0.007D0 !0.007D0 ! shape parameter REAL*8 :: S0, S1 INTEGER, PARAMETER :: COMPSEG = 51 ! NUMBER OF INTERVALS IN SIGMA PROFILE (-0.025 TO 0.025) REAL*8, DIMENSION (:,:,:,:), ALLOCATABLE :: DELTAW REAL*8, DIMENSION (:,:,:), ALLOCATABLE :: SIGMA, SEGGAMMAPURE, SEGGAMMAPUREOLD, CONPR REAL*8, DIMENSION (:,:), ALLOCATABLE :: WHB, SEGGAMMA, SEGGAMMAOLD, CONVPURE REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMASOLUTE_NHB, SIGMASOLUTE_OH, SIGMASOLUTE_OT, SIGMASOLUTE_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_OH, SIGMA1SOLVENT_OT, SIGMA1SOLVENT_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: SIGMA2SOLVENT_NHB, SIGMA2SOLVENT_OH, SIGMA2SOLVENT_OT, SIGMA2SOLVENT_HB REAL*8, DIMENSION (:,:), ALLOCATABLE :: MIXPROFILE, COUNTERSOLUTE REAL*8, DIMENSION (:,:), ALLOCATABLE :: COUNTER1SOLVENT, COUNTER2SOLVENT, CONVERG REAL*8, DIMENSION (:), ALLOCATABLE :: X1DATA, X2DATA, X, COUNTER, X3DATA REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMO, ACOSMO, RNORM, QNORM, DELTAHF REAL*8, DIMENSION (:), ALLOCATABLE :: THETA, PHI, LSG, LNGAMMASG, SUMGAMMA, GAMMA, LNGAMMA, TM, SYSTEMP REAL*8, DIMENSION (:), ALLOCATABLE :: VCOSMOSOLUTE, VCOSMO1SOLVENT , VCOSMO2SOLVENT REAL*8, DIMENSION (:), ALLOCATABLE :: ACOSMOSOLUTE, ACOSMO1SOLVENT , ACOSMO2SOLVENT REAL*8 :: FPOL, ALPHA, ALPHAPRIME, EPS, SUMMATION , CES REAL*8 :: RELTOL, CHANGE, X2OLD, X1NEW, X2NEW, MFSUM, X3OLD INTEGER :: I, J, K, L, P, T, R, ITER, NUMPOINTS, ISTAT , MAXITS, COMP, NUMCOMP, H, B INTEGER, DIMENSION (:), ALLOCATABLE :: UNIT1SOLV, UNITSOLU, UNIT2SOLV CHARACTER (40) :: INPUTFILE, OUTPUTFILE CHARACTER (30), DIMENSION (:), ALLOCATABLE :: SOLVENT1FILE, SOLUTEFILE, SOLVENT2FILE CHARACTER (30), DIMENSION (:), ALLOCATABLE :: SOLV1, SOLV2, SOLU EPS = 3.667D0 ! DIELECTRIC CONSTANT, LIN AND SANDLER USE A CONSTANT FPOL WHICH YIELDS EPS=3.68 !FPOL = (EPS-1.0)/(EPS+0.5) !UNITLESS FPOL = 0.6916D0 !UNITLESS ALPHA = (0.3*AEFFPRIME**(1.5))/(EO) ! KCAL*ANG^4/MOL/e^2 ALPHAPRIME = FPOL*ALPHA ! MISFIT ENERGY CONSTANT, KCAL*ANG^4/MOL/e^2 COMP = 3 !PURE SOLVENT [USER SPECIFIED] ! ALLOCATE INPUT ARRAYS INPUTFILE = "data_10refined_ternaire.txt" ![USER SPECIFIED FILE LOCATION] !NUMPOINTS = 41![USER SPECIFIED, MAXIMUM VALUE: 10,000] !NUMCOMP = 11 ! Pour rentrer manuellement le nombre de points de température ! lecture du nombre de points de température OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD",ACTION="READ", POSITION="REWIND") READ (12,*) NUMCOMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 CLOSE(12) write(*,*) NUMCOMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 ALLOCATE (SOLV1(NUMPOINTS), VCOSMO1SOLVENT(NUMPOINTS), & SOLU(NUMPOINTS), VCOSMOSOLUTE(NUMPOINTS), TM(NUMPOINTS), DELTAHF(NUMPOINTS),& SYSTEMP(NUMPOINTS), X2DATA(NUMPOINTS), UNIT1SOLV(NUMPOINTS),SOLUTEFILE(NUMPOINTS),& SOLVENT1FILE(NUMPOINTS), UNITSOLU(NUMPOINTS), ACOSMO1SOLVENT(NUMPOINTS),& SIGMASOLUTE_NHB(COMPSEG,NUMPOINTS), SIGMASOLUTE_OH(COMPSEG,NUMPOINTS),& SIGMASOLUTE_OT(COMPSEG,NUMPOINTS), SIGMASOLUTE_HB(COMPSEG,NUMPOINTS),& COUNTERSOLUTE(COMPSEG,NUMPOINTS), ACOSMOSOLUTE(NUMPOINTS), & SIGMA1SOLVENT_NHB(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_OH(COMPSEG,NUMPOINTS), & SIGMA1SOLVENT_OT(COMPSEG,NUMPOINTS), SIGMA1SOLVENT_HB(COMPSEG,NUMPOINTS), & COUNTER1SOLVENT(COMPSEG,NUMPOINTS), & ! pour ternaire SIGMA2SOLVENT_NHB(COMPSEG,NUMPOINTS), SIGMA2SOLVENT_HB(COMPSEG,NUMPOINTS), & SIGMA2SOLVENT_OH(COMPSEG,NUMPOINTS), SIGMA2SOLVENT_OT(COMPSEG,NUMPOINTS), & COUNTER2SOLVENT(COMPSEG,NUMPOINTS), & ACOSMO2SOLVENT(NUMPOINTS), SOLV2(NUMPOINTS), VCOSMO2SOLVENT(NUMPOINTS), & UNIT2SOLV(NUMPOINTS),X1DATA(NUMPOINTS), X3DATA(NUMPOINTS), SOLVENT2FILE(NUMPOINTS), & STAT = ISTAT) WRITE (*,*) ISTAT, " ALLOCATION SUCCESSFUL - INPUT ARRAYS" ![USER SPECIFIED INPUT FILE] USE TEMPLATE, L IS COUNTER FOR NUMBER OF POINTS !COLUMN 1: SOLUTE FILE NAME , SOLU(L) !COLUMN 2: SOLUTE COSMO CALCULATION VOLUME, VCOSMOSOLUTE(L) !COLUMN 3: SOLVENT FILE NAME, SOLV1(L) !COLUMN 4: SOLVENT COSMO CALCULATION VOLUME, VCOSMO1SOLVENT(L) !COLUMN 5: MELTING TEMPERATURE, TM(L) !COLUMN 6: LATENT HEAT OF FUSION, DELTAHF(L) en J/mol il me semble !COLUMN 7: SYSTEM TEMPERATURE, SYSTEM(L) !COLUMN 8: SOLUTE MOLE FRACTION, X2DATA(L) ! ######## MANIERE NORMALE DE FAIRE : ############# ! OPEN (UNIT=12, FILE=INPUTFILE, STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(12,*) NUMCOMP, NUMPOINTS, COH_OH, COT_OT, COH_OT, S0, S1 DO L=1, NUMPOINTS READ (12,*) SOLU(L), VCOSMOSOLUTE(L), SOLV1(L), VCOSMO1SOLVENT(L), SOLV2(L), & VCOSMO2SOLVENT(L), TM(L), DELTAHF(L), SYSTEMP(L), X2DATA(L) END DO CLOSE (12) DO I = 1, NUMPOINTS UNIT1SOLV (I) = I + 20 UNIT2SOLV (I) = I + 1020 UNITSOLU (I) = I + 10020 SOLVENT1FILE (I) = TRIM(SOLV1(I))//'_10refined.txt' ![USER SPECIFIED FILE LOCATION] SOLVENT2FILE (I) = TRIM(SOLV2(I))//'_10refined.txt' ![USER SPECIFIED FILE LOCATION] SOLUTEFILE (I) = TRIM(SOLU(I))//'_10refined.txt' ![USER SPECIFIED FILE LOCATION] WRITE (*,*) SOLVENT1FILE(I), SOLVENT2FILE(I), SOLUTEFILE(I), X2DATA(I), TM(I), DELTAHF(I), SYSTEMP(I) END DO ! SIGMA1SOLVENT_HB = 0.0D0 ! SIGMASOLUTE_HB = 0.0D0 DO K = 1, NUMPOINTS OPEN (UNIT=UNIT1SOLV(K), FILE = SOLVENT1FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO1SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT1SOLV(K),*) COUNTER1SOLVENT(J,K), SIGMA1SOLVENT_NHB(J,K), & SIGMA1SOLVENT_OH(J,K), SIGMA1SOLVENT_OT(J,K) ACOSMO1SOLVENT(K) = ACOSMO1SOLVENT(K) + SIGMA1SOLVENT_NHB(J,K) + & SIGMA1SOLVENT_OH(J,K) + SIGMA1SOLVENT_OT(J,K) SIGMA1SOLVENT_HB(J,K) = SIGMA1SOLVENT_OH(J,K) + SIGMA1SOLVENT_OT(J,K) END DO CLOSE (UNIT1SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNIT2SOLV(K), FILE = SOLVENT2FILE(K), STATUS="OLD", ACTION="READ", & POSITION="REWIND") ACOSMO2SOLVENT(K)=0 ! RAJOUT pour initialiser la somme DO J=1, COMPSEG READ(UNIT2SOLV(K),*) COUNTER2SOLVENT(J,K), SIGMA2SOLVENT_NHB(J,K), & SIGMA2SOLVENT_OH(J,K), SIGMA2SOLVENT_OT(J,K) ACOSMO2SOLVENT(K) = ACOSMO2SOLVENT(K) + SIGMA2SOLVENT_NHB(J,K) + & SIGMA2SOLVENT_OH(J,K) + SIGMA2SOLVENT_OT(J,K) SIGMA2SOLVENT_HB(J,K) = SIGMA2SOLVENT_OH(J,K) + SIGMA2SOLVENT_OT(J,K) END DO CLOSE (UNIT2SOLV(K)) END DO DO K = 1, NUMPOINTS OPEN (UNIT=UNITSOLU(K), FILE = SOLUTEFILE(K), STATUS="OLD", ACTION="READ", POSITION = "REWIND") ACOSMOSOLUTE(K)=0 ! RAJOUT pour initialiser la somme DO J = 1, COMPSEG READ(UNITSOLU(K),*) COUNTERSOLUTE(J,K), SIGMASOLUTE_NHB(J,K),& SIGMASOLUTE_OH(J,K), SIGMASOLUTE_OT(J,K) ACOSMOSOLUTE(K) = ACOSMOSOLUTE(K) + SIGMASOLUTE_NHB(J,K) +& SIGMASOLUTE_OH(J,K) + SIGMASOLUTE_OT(J,K) SIGMASOLUTE_HB(J,K) = SIGMASOLUTE_OH(J,K) + SIGMASOLUTE_OT(J,K) END DO CLOSE (UNITSOLU(K)) END DO ![USER SPECIFIED FILE LOCATION] OUTPUTFILE = "results_mullins_10refined.txt" OPEN (UNIT = 13, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_10refined_fractions.txt" OPEN (UNIT = 15, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "results_mullins_10refined_gammas.txt" OPEN (UNIT = 16, FILE = OUTPUTFILE, STATUS = "REPLACE") OUTPUTFILE = "control.txt" OPEN (UNIT = 17, FILE = OUTPUTFILE, STATUS = "REPLACE") 2 FORMAT (1X,A1,1X,A1, 1X, A16, 1X, A16) 3 FORMAT (1X,A5,1X,A4,1X,A16,1X,A16,1X,A4,1X,A2,1X,A2,1X,A7,1X,A7) !4 format (1F8.4) WRITE (13,2) "-", "-", "COMP(1)", "COMP(2)" WRITE (13,3) "POINT","ITER","SOLU","SOLV1","TEMP","X1","X2","X3" !ALLOCATE ARRAYS USED IN MAIN ITERATIVE LOOP ALLOCATE (SIGMA(COMPSEG,COMP, 3), COUNTER(COMPSEG),& MIXPROFILE(COMPSEG,3), DELTAW(COMPSEG,COMPSEG,3,3),& SEGGAMMA(COMPSEG,3),SEGGAMMAOLD(COMPSEG,3),CONVERG(COMPSEG,3),SEGGAMMAPURE(COMPSEG,COMP,3),& SEGGAMMAPUREOLD(COMPSEG,COMP,3), CONVPURE(COMPSEG,COMP), WHB(COMPSEG, COMPSEG), & VCOSMO(COMP), ACOSMO(COMP), X(COMP), THETA(COMP),PHI(COMP),LNGAMMASG(COMP),GAMMA(COMP), & LSG(COMP), SUMGAMMA(COMP), LNGAMMA(COMP),RNORM(COMP), QNORM(COMP), & CONPR(COMPSEG, COMP, 3), STAT = ISTAT) WRITE(*,*) ISTAT, " ALLOCATION SUCCESSFUL - CALCULATION ARRAYS" ! Begin iteration loop for all binaries DO P=1, NUMPOINTS ! write(*,*) SYSTEMP(P) ! BEGIN ITERATION LOOP FOR ALL POINTS (composition), REASSIGNING NEW SIGMA PROFILES FOR EACH ITERATIONS (pas besoin) DO R = 1, NUMCOMP ! RE-ZERO CALCULATION LOOP ARRAYS SIGMA = 0.0D0 COUNTER = 0.0D0 MIXPROFILE = 0.0D0 DELTAW = 0.0D0 SEGGAMMA = 1.0D0 SEGGAMMAOLD = 1.0D0 CONVERG = 1.0D0 SEGGAMMAPURE = 1.0D0 SEGGAMMAPUREOLD = 1.0D0 WHB = 0.0D0 THETA = 0.0D0 PHI = 0.0D0 LSG = 0.0D0 RNORM = 0.0D0 QNORM = 0.0D0 LNGAMMASG = 0.0D0 SUMGAMMA = 0.0D0 LNGAMMA = 0.0D0 GAMMA = 0.0D0 ACOSMO = 0.0D0 VCOSMO = 0.0D0 ! EXPRESSION OF CES CES = AES + BES/(SYSTEMP(P)**2) !CES = ALPHAPRIME/2 DO I = 1, COMPSEG !COUNTER(I) = COUNTER1SOLVENT(I,P) COUNTER(I) = COUNTER1SOLVENT(I,P) ! 1->NHB, 2->OH, 3->OT ! SIGMA(SEGMENT, SPECIES, NHB or OH or OT) SIGMA(I,1,1) = SIGMA1SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMA1SOLVENT_NHB(I,P) SIGMA(I,2,1) = SIGMA2SOLVENT_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S0**2)))& + SIGMA2SOLVENT_NHB(I,P) SIGMA(I,3,1) = SIGMASOLUTE_HB(I,P)*(dexp(-COUNTER(I)**2/(2*S1**2)))& + SIGMASOLUTE_NHB(I,P) SIGMA(I,1,2) = SIGMA1SOLVENT_OH(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,2,2) = SIGMA2SOLVENT_OH(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,3,2) = SIGMASOLUTE_OH(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S1**2))) SIGMA(I,1,3) = SIGMA1SOLVENT_OT(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,2,3) = SIGMA2SOLVENT_OT(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S0**2))) SIGMA(I,3,3) = SIGMASOLUTE_OT(I,P)*(1.0D0-dexp(-COUNTER(I)**2/(2*S1**2))) ! SIGMA(I,1,1) = SIGMA1SOLVENT_NHB(I,P) ! SIGMA(I,2,1) = SIGMASOLUTE_NHB(I,P) ! SIGMA(I,1,2) = SIGMA1SOLVENT_OH(I,P) ! SIGMA(I,2,2) = SIGMASOLUTE_OH(I,P) ! SIGMA(I,1,3) = SIGMA1SOLVENT_OT(I,P) ! SIGMA(I,2,3) = SIGMASOLUTE_OT(I,P) END DO ACOSMO(1) = ACOSMO1SOLVENT(P) ACOSMO(2) = ACOSMO2SOLVENT(P) ACOSMO(3) = ACOSMOSOLUTE(P) VCOSMO(1) = VCOSMO1SOLVENT(P) VCOSMO(2) = VCOSMO2SOLVENT(P) VCOSMO(3) = VCOSMOSOLUTE(P) X(1) = 1.0D0 - X2DATA(P) X(2) = X2DATA(P) X(3) = 0.000000001D0 !write(*,*) ACOSMO(1), ACOSMO(2), VCOSMO(1), VCOSMO(2), X(1), X(2) ! CONVERGENCE LOOP, MAX ITERATIONS = 1000 DO ITER = 1, 1000 X(1) = 1.0D0 - X(2) - X(3) X2OLD = X(2) X3OLD = X(3) ! CALCULATE SIGMA PROFILE FOR MIXTURE FOR EACH DATA POINT ! Ps(SIGMA) = SUM(Xi*Ai*Pi(SIGMA))/SUM(Xi*Ai) ! SIGMA(J,K) = P'(SIGMA) = Pi(SIGMA)*Ai !WRITE (13,*) " MIXPROFILE" DO J = 1,COMPSEG ! 51 SEGMENTS IN MIXTURE SIGMA PROFILE ![MIXTURE PROFILE FOR PURE SOLVENT SYSTEM] ! MIXPROFILE IS p(sigmai(J)) FOR THE MIXTURE in case NHB (1), OH (2), OT(3) MIXPROFILE(J,1) = ((X(1)*SIGMA(J,1,1))+(X(2)*SIGMA(J,2,1)) + (X(3)*SIGMA(J,3,1))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) MIXPROFILE(J,2) = ((X(1)*SIGMA(J,1,2))+(X(2)*SIGMA(J,2,2))+(X(3)*SIGMA(J,3,2))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) MIXPROFILE(J,3) = ((X(1)*SIGMA(J,1,3))+(X(2)*SIGMA(J,2,3))+(X(3)*SIGMA(J,3,3))) / & ((X(1)*ACOSMO(1)) + (X(2)*ACOSMO(2)) + (X(3)*ACOSMO(3))) END DO ! DETERMINE THE LIN AND SANDLER EXCHANGE ENERGY, DELTAW, KCAL/MOL (ORIGINAL MODEL) DO H=1,3 DO B=1,3 DO I = 1, COMPSEG DO K = 1, COMPSEG IF (H==B.AND.H==2.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN ! DELTAW(I,K,H,B) is delta for segment I for H=1 (NHB) or H=2 (OH) H=3 (OT) for segment K DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COH_OH*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (H==B.AND.H==3.AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COT_OT*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (((H==2.AND.B==3).OR.(H==3.AND.B==2)).AND.(COUNTER(I)*COUNTER(K))<0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 - & COH_OT*(COUNTER(I)-COUNTER(K))**2 ENDIF IF (H==1.OR.B==1.OR.(COUNTER(I)*COUNTER(K))>=0.0D0) THEN DELTAW(I,K,H,B) = CES*(COUNTER(I)+COUNTER(K))**2.0D0 ENDIF END DO END DO ENDDO ENDDO ! ITERATION FOR MIXTURE SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 SEGGAMMA = 1.0D0 !INTIAL VALUE FOR MIXTURE SEGMENT ACTIVITY COEF. DO SEGGAMMAOLD = SEGGAMMA DO B=1,3 DO I = 1, COMPSEG SUMMATION = 0.0D0 DO H=1,3 DO K = 1, COMPSEG SUMMATION = SUMMATION + MIXPROFILE(K,H)*SEGGAMMAOLD(K,H) * & DEXP(-DELTAW(I,K,B,H)/(RGAS * SYSTEMP(P))) ENDDO ENDDO SEGGAMMA(I,B) = DEXP(-DLOG(SUMMATION)) SEGGAMMA(I,B) = (SEGGAMMA(I,B) + SEGGAMMAOLD(I,B))/2.0D0 END DO ENDDO DO B=1, 3 DO I=1, COMPSEG CONVERG(I,B) = DABS((SEGGAMMA(I,B) - SEGGAMMAOLD(I,B))/SEGGAMMAOLD(I,B)) END DO ENDDO IF (MAXVAL(CONVERG) <=0.000001) EXIT ENDDO ! ITERATION FOR PURE SPECIES, L, SEGMENT ACTIVITY COEF. ! CONVERGENCE BY DAMPING METHOD, DAMPING FACTOR = 1/2 ! ITERATION FOR SEGMENT ACITIVITY COEF (PURE SPECIES) ! HERE p(sigma) is just SIGMA/Ai total IF (ITER==1) THEN DO L = 1, COMP !SEGGAMMAPURE (:,L,:) = 1.0D0 SEGGAMMAPURE = 1.0D0 DO !SEGGAMMAPUREOLD (:,L,:) = SEGGAMMAPURE (:,L,:) SEGGAMMAPUREOLD = SEGGAMMAPURE DO I = 1, COMPSEG DO B=1,3 SUMMATION = 0.0D0 DO H=1,3 DO K = 1, COMPSEG SUMMATION = SUMMATION + (SIGMA(K,L,H)/ACOSMO(L)) * & SEGGAMMAPUREOLD(K,L,H)*DEXP(-DELTAW(I,K,B,H)/(RGAS*SYSTEMP(P))) ENDDO END DO SEGGAMMAPURE(I,L,B)=DEXP(-DLOG(SUMMATION)) SEGGAMMAPURE(I,L,B)=(SEGGAMMAPURE(I,L,B)+SEGGAMMAPUREOLD(I,L,B))/2.0D0 ENDDO END DO DO B=1,3 DO I=1,COMPSEG CONPR(I,L,B)=DABS((SEGGAMMAPURE(I,L,B)-SEGGAMMAPUREOLD(I,L,B))/& SEGGAMMAPUREOLD(I,L,B)) ENDDO ENDDO IF (MAXVAL(CONPR) <=0.000001) EXIT ENDDO ENDDO ENDIF ! THE STAVERMAN-GUGGENHEIM EQUATION DO I = 1,COMP RNORM(I) = VCOSMO(I)/VNORM QNORM(I) = ACOSMO(I)/ANORM END DO DO I = 1, COMP THETA(I) = (X(I)*QNORM(I))/(X(1)*QNORM(1) + X(2)*QNORM(2) + X(3)*QNORM(3)) PHI(I) = (X(I)*RNORM(I))/(X(1)*RNORM(1) + X(2)*RNORM(2) + X(3)*RNORM(3)) LSG(I) = (COORD/2.0D0)*(RNORM(I)-QNORM(I))-(RNORM(I)-1.0D0) END DO ! GAMMASG1 AND GAMMASG2 ARE ACTUALLY LNGAMMASG DO I = 1, COMP LNGAMMASG(I) = DLOG(PHI(I)/X(I)) + (COORD/2.0D0)*QNORM(I)* & DLOG(THETA(I)/PHI(I)) + LSG(I) - (PHI(I)/X(I))*(X(1)*LSG(1) + X(2)*LSG(2) + X(3)*LSG(3)) END DO !CALCULATION OF GAMMAS SUMGAMMA = 0.0D0 DO K = 1, COMP DO B = 1, 3 DO I=1, COMPSEG SUMGAMMA(K) = SUMGAMMA(K)+((SIGMA(I,K,B)/AEFFPRIME)*(DLOG(SEGGAMMA(I,B)/ & SEGGAMMAPURE(I,K,B)))) END DO ENDDO GAMMA(K) = DEXP(SUMGAMMA(K) + (LNGAMMASG(K))) LNGAMMA(K) = DLOG(GAMMA(K)) END DO ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(3) = DEXP((DELTAHF(P)/8.3145D0/TM(P))*(1.0D0-TM(P)/SYSTEMP(P)) - LNGAMMA(3)) ! Affichage des résultats intermédiaires ! WRITE(*,*) P, X(2) !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(3) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLE FRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(3) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT ! DAMPING FACTOR (1/3) USED IN CONVERGING X(2) X(3) = (X(3) + 2.0D0*X3OLD)/3.0D0 MFSUM = X(1)+X(2)+X(3) X(1) = X(1)/MFSUM X(2) = X(2)/MFSUM X(3) = X(3)/MFSUM ! Pour afficher la fraction molaire finale intermédiaire ! write(*,*) 'final', X(2) ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = DABS((X(3) - X3OLD)/(X3OLD + 1.0D-16)) !WRITE (*,*) " RELATIVE CHANGE:", P, CHANGE ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) GOTO 10 write(17,*) X3OLD, X(3), LNGAMMASG(3), SUMGAMMA(1) !END ITERATIVE LOOP END DO ! POINT, ITER, SOLU, SOLV1,SYSTEMP, X1, X2, GAMMA1, GAMMA2 (INPUT FORMAT TXT FILE) !4 FORMAT (1X,I4,2X,I4,2X,A16,2X,A16,2X,E,1X,E,1X,E,1X,E,1X,E) 10 WRITE (13,*) P,ITER,SOLU(P),SOLV1(P),SYSTEMP(P),X(1),X(2),GAMMA(1),GAMMA(2) ! END OF CALCULATIONS FOR ALL POINTS write (15,*) X(3) write (16,*) GAMMA(3) write(*,*) X(3) ! #### CAS OU ON MODIFIE LA COMPOSITION : write(*,*) "pour X2DATA=", X2DATA(P) X2DATA(P) = X2DATA(P) - 0.10 END DO ! Iteration composition (variable R) END DO ! Iteration points (variable P) CLOSE (13) ! CLOSE OUTPUT FILE close (15) ! CLOSE l'AUTRE OUTPUT FILE close (16) ! CLOSE l'AUTRE OUTPUT FILE ! DEALLOCATE ARRAYS FOR NEXT POINT DEALLOCATE (SIGMA, COUNTER, MIXPROFILE, DELTAW, & SEGGAMMA, SEGGAMMAOLD, CONVERG, SEGGAMMAPURE, & SEGGAMMAPUREOLD, WHB, CONVPURE, & VCOSMO, ACOSMO, X, THETA, PHI, LNGAMMASG, GAMMA, & LSG, SUMGAMMA, LNGAMMA, RNORM, QNORM, CONPR, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - CALCULATION ARRAYS" DEALLOCATE (SOLV1, VCOSMO1SOLVENT, SOLU, VCOSMOSOLUTE, TM, DELTAHF, SYSTEMP, & X2DATA, UNIT1SOLV, SOLUTEFILE, SOLVENT1FILE, UNITSOLU, & SIGMASOLUTE_NHB, SIGMASOLUTE_OH, SIGMASOLUTE_OT, SIGMASOLUTE_HB, & COUNTERSOLUTE, ACOSMOSOLUTE, SIGMA1SOLVENT_NHB, SIGMA1SOLVENT_OH, & COUNTER1SOLVENT, SIGMA1SOLVENT_OT, SIGMA1SOLVENT_HB, ACOSMO1SOLVENT, & ACOSMO2SOLVENT, VCOSMO2SOLVENT, UNIT2SOLV, SOLVENT2FILE, SOLV2, & SIGMA2SOLVENT_NHB, SIGMA2SOLVENT_OH, SIGMA2SOLVENT_OT, X3DATA, X1DATA, & SIGMA2SOLVENT_HB, COUNTER2SOLVENT, STAT = ISTAT) WRITE (*,*) ISTAT, " DEALLOCATION SUCCESSFUL - INPUT ARRAYS" END PROGRAM GAMMASOLUBILITY II/NRTL-SAC II.1/NRTL-SAC binaire II.1.1/NRTL-SAC binaire FORTRAN PROGRAM NRTL_SAC !****************************************************************************** ! GENERAL COMMENTS: ! Author: Adel Ghaderi ! Supervisor: Eric Mullins ! Advisor: Y. A. Liu ! Institution: Virginia Polytechnic Institute and State University ! Date: October 27,2006 ! Code is setup to handle pure solvent systems for solubility of single solute ! for subscripts referring to the segments ! 1 refers to X ! 2 refers to Y- ! 3 refers to Y+ ! 4 refers to Z ! Program is coded for 4 segments; simply changing value of NUMSEG is not ! sufficient to vary number of segments. !****************************************************************************** ! REFERENCES: ! Chen, C.C.; Song, Y.; "Solubility modeling with a nonrandom two-liquid segment activity coefficient ! model," IECR (2004), v.43, no.26, 8354-8362. !****************************************************************************** IMPLICIT NONE INTEGER, PARAMETER :: SINGLE = 4 INTEGER, PARAMETER :: DOUBLE = 8 INTEGER :: i, j, k, m, n, y, z, ITER, ITERcount, NUMSPEC, NUMPOINTS, NUMSYS, istat INTEGER, PARAMETER :: NUMSEG = 4 !number of segments; 4 (X, Y-, Y+, Z) as DEFINED by C. C. Chen INTEGER, PARAMETER :: CHARLEN = 100 !maximum allowable length name of the system REAL(KIND=DOUBLE) :: LNSACTERM1TOP, LNSACTERM1BOT, LNSACTERM2 REAL(KIND=DOUBLE) :: LNSACITERM1TOP, LNSACITERM1BOT, LNSACITERM2 REAL(KIND=DOUBLE) :: LNSACTERM2BOT, LNSACTERM4TOP, LNSACTERM4BOT REAL(KIND=DOUBLE) :: LNSACITERM2BOT, LNSACITERM4TOP, LNSACITERM4BOT REAL(KIND=DOUBLE) :: LNSAC, LNSACI, dummy, XOLD, MFSUM, RELTOL, CHANGE, LNRE, SOM REAL(KIND=DOUBLE), ALLOCATABLE, DIMENSION(:) :: T, XSEG, LNGAMMAR, LNGAMMAC, PHI, dHfus, TM, Rt REAL(KIND=DOUBLE), ALLOCATABLE, DIMENSION(:,:) :: TAU, ALPHA, G, R, XORIG, X, XSPEC, GAMMA CHARACTER(CHARLEN), ALLOCATABLE, DIMENSION(:) :: TITLE CHARACTER(CHARLEN), ALLOCATABLE, DIMENSION(:,:) :: NAME !------------------------------------------------------------------------------------ ! Initialize Tau and Alpha values, from C.C. Chen, et al. !------------------------------------------------------------------------------------ ALLOCATE( TAU(NUMSEG,NUMSEG), ALPHA(NUMSEG,NUMSEG), G(NUMSEG,NUMSEG) ) DO i = 1, 4 ! 4 segment species (X, Y-, Y+, Z) DO j = 1, 4 TAU(i,j) = 0 ALPHA(i,j) = 0 ! note ALPHA(i,j) = ALPHA(j,i) END DO END DO TAU(1,2) = 1.643 ! 1 denotes X (hydrophobic segment) TAU(2,1) = 1.834 ! 2 denotes Y- (polar attractive segment) ALPHA(1,2) =0.200; ALPHA(2,1) = .200 TAU(1,4) = 6.547 ! 4 denotes Z (hydrophilic segment) TAU(4,1) = 10.949 ALPHA(1,4) =0.200; ALPHA(4,1) = .200 TAU(2,4) = -2.000 TAU(4,2) = 1.787 ALPHA(2,4) =0.300; ALPHA(4,2) = .300 TAU(3,4) = 2.000 ! 3 denotes Y+ (polar repulsive segment) TAU(4,3) = 1.787 ALPHA(3,4) =0.300; ALPHA(4,3) = .300 TAU(1,3) = 1.643 TAU(3,1) = 1.834 ALPHA(1,3) =0.200; ALPHA(3,1) = .200 !------------------------------------------------------------------------------------ ! Calculate values for G DO i = 1, 4 DO j = 1, 4 G(i,j) = exp( - ALPHA(i,j) * TAU(i,j) ) !write(*,*) 'G(',i,',',j,') =',G(i,j) END DO END DO ! Generate input file OPEN(UNIT = 13, FILE = "base_input.txt", & STATUS="OLD", ACTION="READ", POSITION="REWIND") ! File will be closed at end of section OPEN(UNIT=14,FILE="input_intermediaire.txt",STATUS="REPLACE", ACTION="WRITE") ! File will be closed at end of section ! The input data of interest is for pure component solubility NUMSPEC = 2 ! For simplicity, every point is treated as a single system (this is usually the case anyway) NUMPOINTS = 1 119 FORMAT( I4 ) 120 FORMAT( I1, 5X, I1 ) 121 FORMAT( A15 ) 122 FORMAT( I1,2X,F10.2,2X,F7.2,2X,F12.10,2X,F12.10,2X,F12.10,2X,F12.10 ) 123 FORMAT( F8.2, 2X, F17.15, 2X, F17.15 ) 124 FORMAT( I1,2X,I1,2X,I1,2X,F12.10,2X,F12.10,2X,F12.10,2X,F12.10 ) READ(13,119) NUMSYS !NUMSYS = 1 WRITE(14,119) NUMSYS ALLOCATE( R(8,NUMSYS), dHfus(NUMSYS), TM(NUMSYS), T(NUMSYS), X(NUMSYS,2), NAME(NUMSYS,2), STAT = istat ) WRITE(*,*) istat , " ALLOCATION SUCCESSFUL IF 0" DO i = 1, NUMSYS ! INPUT FILE FORMAT: SOLUTE NAME, HEAT OF FUSION, MELTING TEMPERATURE, SEGMENT PARAMETERS (X,Y-, ! Y+,Z),SOLVENT NAME, SEGMENT PARAMETERS (X,Y-,Y+,Z), SYSTEM TEMPERATURE, EXPERIMENTAL SOLUBILITY READ(13,*) NAME(i,2), dHfus(i), TM(i), R(1,i), R(2,i), R(3,i), R(4,i), NAME(i,1), R(5,i), R(6,i), & &R(7,i), R(8,i), T(i), X(i,2) X(i,1) = 1.0D0 - X(i,2) WRITE(14,120) NUMSPEC, NUMPOINTS WRITE(14,121) NAME(i,1) WRITE(14,121) NAME(i,2) WRITE(14,124) 1, 0, 0, R(5,i), R(6,i), R(7,i), R(8,i) WRITE(14,122) 2, dHfus(i), TM(i), R(1,i), R(2,i), R(3,i), R(4,i) WRITE(14,123) T(i), X(i,1), X(i,2) END DO DEALLOCATE( R, dHfus, TM, T, X, NAME, STAT = istat ) WRITE(*,*) istat, " DEALLOCATION SUCCESSFUL IF 0" CLOSE(13) CLOSE(14) !------------------------------------------------------------------------------------ ! Open input and output files OPEN(UNIT = 11, FILE = "input_intermediaire.txt",STATUS="OLD", ACTION="READ", POSITION="REWIND") ! File will be closed at end of program READ(11,*) NUMSYS OPEN(UNIT=12, FILE = "output.txt",& STATUS="REPLACE", ACTION="WRITE") OPEN(UNIT=51, FILE = "solubilite.txt",& STATUS="REPLACE", ACTION="WRITE") ! File will be closed at end of program WRITE(12,115) "ITER", "SOLUTE", "SOLVENT", "TEMP(K)", "GAMMA-1", "GAMMA-2", "X1-EXP", "X2-exp", "X1-pred",& "X2-pred", "Rel.-%E-(X2)", "%LN(E)" WRITE(12,115) "====", "======", "=======", "=======", "=======", "=======", "======", "======", "=======",& "=======", "============", "======" 115 FORMAT(2X, A4, 3X, A6, 6X, A7, 4X, A7, 4X, A7, 8X, A7, 8X, A6, 6X, A6, 7X, A7, 6X, A7, 6X,& &A12, 5X, A6) ! Iterate over all different multicomponent systems !------------------------------------------------------------------------------------ DO z = 1, NUMSYS ! Some file manipulation and memory allocation READ(11,*) NUMSPEC, NUMPOINTS !number of system components and number of data points; used in memory allocation WRITE(*,*) NUMSPEC, NUMPOINTS ALLOCATE( NAME(NUMSPEC,1), TITLE(NUMSYS), dHfus(NUMSPEC), TM(NUMSPEC), & T(NUMPOINTS), XSEG(NUMSEG), LNGAMMAR(NUMSPEC), LNGAMMAC(NUMSPEC), PHI(NUMSPEC), & R(NUMSEG,NUMSPEC), XSPEC(NUMSEG,NUMSPEC), XORIG(NUMPOINTS,NUMSPEC), X(NUMPOINTS,NUMSPEC), & GAMMA(NUMPOINTS,NUMSPEC), Rt(NUMSPEC), STAT = istat ) WRITE(*,*) istat , " ALLOCATION SUCCESSFUL IF 0" ! Collect VLE data from file !READ(11,*) TITLE(z) DO i = 1, NUMSPEC READ(11,*) NAME(i,1) END DO DO i = 1, NUMSPEC READ(11,*) j , dHfus(i), TM(i), R(1,i), R(2,i), R(3,i), R(4,i) ! values for X, Y-, Y+, and Z for species i END DO IF ( NUMSPEC == 2 ) THEN DO i = 1, NUMPOINTS READ(11,*) T(i), XORIG(i,1), XORIG(i,2) ! PURE Solvent X(i,1) = XORIG(i,1) X(i,2) = XORIG(i,2) END DO ELSE IF ( NUMSPEC == 3 ) THEN DO i = 1, NUMPOINTS READ(11,*) T(i), XORIG(i,1), XORIG(i,2), XORIG(i,3) ! BINARY Solvent X(i,1) = XORIG(i,1) X(i,2) = XORIG(i,2) X(i,3) = XORIG(i,3) END DO ELSE IF ( NUMSPEC == 4 ) THEN DO i = 1, NUMPOINTS READ(11,*) T(i), XORIG(i,1), XORIG(i,2), XORIG(i,3), XORIG(i,4) ! TERNARY Solvent X(i,1) = XORIG(i,1) X(i,2) = XORIG(i,2) X(i,3) = XORIG(i,3) X(i,4) = XORIG(i,4) END DO ELSE IF ( NUMSPEC == 5 ) THEN DO i = 1, NUMPOINTS !QUATERNARY Solvent READ(11,*) T(i), XORIG(i,1), XORIG(i,2), XORIG(i,3), XORIG(i,4), XORIG(i,5) X(i,1) = XORIG(i,1) X(i,2) = XORIG(i,2) X(i,3) = XORIG(i,3) X(i,4) = XORIG(i,4) X(i,5) = XORIG(i,5) END DO ELSE WRITE(*,*) "The number of components is not 2, 3, 4, or 5; please refer to source code." WRITE(*,*) NUMSPEC PAUSE END IF ! Generate XSPEC, denoted Xj,I in C. C. Chen (note: XSPEC is independent of system composition) et somme R(i,j) !------------------------------------------------------------------------------------ DO i = 1, NUMSPEC dummy = 0 DO j = 1, 4 XSPEC(j,i) = 0 Rt(i)=0 DO k = 1, 4 Rt(i) = Rt(i) + R(k,i) END DO XSPEC(j,i) = R(j,i) / Rt(i) END DO END DO ! calculate over all points DO n = 1, NUMPOINTS DO ITER = 1, 200 ! Generate XSEG, denoted Xj in C. C. Chen !------------------------------------------------------------------------------------ DO i = 1, 4 dummy = 0 XSEG(i) = 0 DO j = 1, NUMSPEC dummy = dummy + X(n,j)*Rt(j) XSEG(i) = XSEG(i) + X(n,j)*R(i,j) END DO XSEG(i) = XSEG(i) / dummy END DO ! calculate over all species DO i = 1, NUMSPEC ! Calculate residual term, LNGAMMAR, from NRTL-SAC model !------------------------------------------------------------------------------------ LNGAMMAR(i) = 0 DO m = 1, 4 LNSACTERM1TOP = 0 LNSACTERM1BOT = 0 LNSACTERM2 = 0 LNSACITERM1TOP = 0 LNSACITERM1BOT = 0 LNSACITERM2 = 0 DO j = 1, 4 !j is same as j, k, and m' for terms 1 and 2 LNSACTERM2BOT = 0 LNSACTERM4TOP = 0 LNSACTERM4BOT = 0 LNSACITERM2BOT = 0 LNSACITERM4TOP = 0 LNSACITERM4BOT = 0 DO k = 1, 4 !k is same as k, j, and k for terms 2 and 4 LNSACTERM2BOT = LNSACTERM2BOT + XSEG(k)*G(k,j) LNSACTERM4TOP = LNSACTERM4TOP + XSEG(k)*G(k,j)*TAU(k,j) LNSACTERM4BOT = LNSACTERM4BOT + XSEG(k)*G(k,j) ! note same as LNSACTERM2BOT LNSACITERM2BOT = LNSACITERM2BOT + XSPEC(k,i)*G(k,j) LNSACITERM4TOP = LNSACITERM4TOP + XSPEC(k,i)*G(k,j)*TAU(k,j) LNSACITERM4BOT = LNSACITERM4BOT + XSPEC(k,i)*G(k,j) ! note same as LNSACITERM2BOT END DO ! end iteration of k !note that G(m,j) in LNSACTERM2 and LNSACITERM2 have indicies reversed !from those of other G values, but this is confirmed by lit (C.C. Chen) LNSACTERM1TOP = LNSACTERM1TOP + XSEG(j)*G(j,m)*TAU(j,m) LNSACTERM1BOT = LNSACTERM1BOT + XSEG(j)*G(j,m) LNSACTERM2 = LNSACTERM2 + XSEG(j) * G(m,j) / & LNSACTERM2BOT * ( TAU(m,j) - LNSACTERM4TOP / LNSACTERM4BOT ) LNSACITERM1TOP = LNSACITERM1TOP + XSPEC(j,i)*G(j,m)*TAU(j,m) LNSACITERM1BOT = LNSACITERM1BOT + XSPEC(j,i)*G(j,m) LNSACITERM2 = LNSACITERM2 + XSPEC(j,i) * G(m,j) / & LNSACITERM2BOT * ( TAU(m,j) - LNSACITERM4TOP / LNSACITERM4BOT ) END DO ! end iteration of j LNSAC = LNSACTERM1TOP / LNSACTERM1BOT + LNSACTERM2 LNSACI = LNSACITERM1TOP / LNSACITERM1BOT + LNSACITERM2 !if (i==2) then ! write(*,*) 'm=',m,'lnac =', LNSAC, 'lnaci =', LNSACI !end if LNGAMMAR(i) = LNGAMMAR(i) + R(m,i) * ( LNSAC - LNSACI ) END DO ! end iteration of m !if (i==2) then !write(*,*) 'ln gamma r�sid =', LNGAMMAR(i) !end if !------------------------------------------------------------------------------------ ! Calculate combinatorial term, LNGAMMAC, from Flory-Huggins, as done by C. C. Chen !------------------------------------------------------------------------------------ ! calcul de PHI(k) DO k = 1, NUMSPEC PHI(k) = 0 DO j = 1, NUMSPEC PHI(k) = PHI(k) + Rt(j)*X(n,j) END DO PHI(k) = Rt(k)*X(n,k)/PHI(K) !write(*,*) Rt(1), Rt(2), 'PHI(',k,') =', PHI(k) END DO ! calcul de SOM = somme sur k de [PHI(k)/Rt(k)] SOM = 0 DO j = 1, NUMSPEC SOM = SOM + PHI(j)/Rt(j) END DO ! calcul de la partie combinatoire, model as shown by C. C. Chen LNGAMMAC(i) = LOG( PHI(i) / X(n,i) ) + 1 - Rt(i) * SOM !write(*,*) 'SOM =', SOM, 'ln gamma r�sid', i, '=', LNGAMMAC(i) ! model as found in G. M. Kontogeorgis et al., Fluid Phase Equilibria 92 (1994) 35-66 ! LNGAMMAC(i) = LOG( PHI(i) / X(n,i) ) + 1 - ( PHI(i) / X(n,i) ) ! note LOG denotes natural log in FORTRAN !------------------------------------------------------------------------------------ ! Calculate overall activity coefficient, GAMMA, and predicted Y, YPRED !------------------------------------------------------------------------------------ GAMMA(n,i) = EXP( LNGAMMAR(i) + LNGAMMAC(i) ) !------------------------------------------------------------------------------------ END DO ! end iteration of i ! write(25,*) 'lngammac =', LNGAMMAC(2) ! XOLD = X(n,NUMSPEC) ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(n,NUMSPEC)=EXP(dHfus(NUMSPEC)/(TM(NUMSPEC)*8.314)*(1-TM(NUMSPEC)/T(n))- & LOG(GAMMA(n,NUMSPEC))) !X(n,NUMSPEC)=EXP(-2.630-LOG(GAMMA(n,NUMSPEC))) !-2.630 = lnKsp ! 8.314 J/mol*K; 0.001 kJ/ 1 J conversion because dHfus is in kJ/mol ! DAMP CALCULATED SOLUTE MOLE FRACTION X(n,NUMSPEC) = ( 2.0D0*XOLD + X(n,NUMSPEC) ) / 3.0D0 !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(n,NUMSPEC) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLEFRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(n,NUMSPEC) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT MFSUM = 0 DO y = 1, NUMSPEC MFSUM = MFSUM + X(n,y) END DO DO y = 1, NUMSPEC X(n,y) = X(n,y) / MFSUM END DO ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = ABS( ( X(n,NUMSPEC) - XOLD)/( XOLD + 1.0D-16 ) ) !WRITE (*,*) " RELATIVE CHANGE:", CHANGE WRITE(*,*) ITER ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) THEN ITERcount = ITER GOTO 200 END IF !END ITERATIVE LOOP END DO ! end iteration of ITER 200 dummy = ABS( X(n,NUMSPEC) - XORIG(n,NUMSPEC) ) / XORIG(n,NUMSPEC) * 100 ! relative error IF ( NUMSPEC == 2 ) THEN write(*,*) ITER WRITE(12,111) ITER, NAME(2,1), NAME(1,1), T(n), GAMMA(n,1), GAMMA(n,2), XORIG(n,1),& XORIG(n,NUMSPEC), X(n,1), X(n,NUMSPEC), dummy ! Pure Solvent !test d'�criture dans le shell WRITE(*,111) ITER, NAME(2,1), NAME(1,1), T(n), GAMMA(n,1), GAMMA(n,2), XORIG(n,1),& XORIG(n,NUMSPEC), X(n,1), X(n,NUMSPEC), dummy ! Pure Solvent !write(*,*) 'lngammar =',LNGAMMAR(2) !write(*,*) 'lngammac =', LNGAMMAC(2) !write(*,*) T(n) write(unit=51,fmt='(F7.5)') X(n,2) 111 FORMAT(2X,I3,4X,A10,1X,A9,1X,F8.2,2X,F12.8,2X,F12.8,5X,F8.5,5X,F8.5,5X,F8.5,5X,F8.5,6X,F8.5) !111 FORMAT(2X,I3,5X,A15,2X,A15,2X,E4.4,2X,E4.4,2X,E4.4,2X,E4.4,2X,E4.4,2X,E4.4,2X,E4.4,2X,E4.4) ELSE WRITE(*,*) "The number of components is not 2 (pure solvent solubility); & &source code must be changed for different systems." PAUSE END IF dummy = 0 END DO ! end iteration of n DEALLOCATE( NAME, TITLE, dHfus, TM, T, XSEG, LNGAMMAR, LNGAMMAC, PHI, R, XSPEC, XORIG, X, GAMMA, STAT& = istat) WRITE(*,*) istat, " DEALLOCATION SUCCESSFUL IF 0" WRITE(*,*) "System ", z, " complete." write(*,*) NUMSYS END DO ! end iteration of z CLOSE(12) CLOSE(11) DEALLOCATE ( TAU, ALPHA, G ) WRITE(*,*) "DONE!" PAUSE END PROGRAM NRTL_SAC II.1.2/NRTL-SAC binaire MATLAB function [y] = nrtlsac1(X) % % nrtlsac1 est une résolution de la méthode NRTLS-SAC % Il renvoie la solubilité prédite (VECTEUR) pour chaque cas considéré % Il lui faut en entrée [X Y- Y+ Z] du soluté et le document input2.txt % qui contient les données expérimentales (une ligne pour chaque) % Ces lignes contiennent : % X Y- Y+ Z (du solvant) Xexp T Tfus deltaHfus % TAU = zeros(4,4); % allocation et initialisation de la matrice TAU ALPHA=zeros(4,4); % allocation et initialisation de la matrice TAU TAU(1,2) = 1.643; % 1 denotes X (hydrophobic segment) TAU(2,1) = 1.834 ; % 2 denotes Y- (polar attractive segment) ALPHA(1,2) =0.200; ALPHA(2,1) = .200; TAU(1,4) = 6.547; % 4 denotes Z (hydrophilic segment) TAU(4,1) = 10.949; ALPHA(1,4) =0.200; ALPHA(4,1) = .200; TAU(2,4) = -2.000; TAU(4,2) = 1.787; ALPHA(2,4) =0.300; ALPHA(4,2) = .300; TAU(3,4) = 2.000; % 3 denotes Y+ (polar repulsive segment) TAU(4,3) = 1.787; ALPHA(3,4) =0.300; ALPHA(4,3) = .300; TAU(1,3) = 1.643; TAU(3,1) = 1.834; ALPHA(1,3) =0.200; ALPHA(3,1) = 0.200; G=zeros(4,4); XSPEC = zeros(4,2); % les Xj,I de la publi avec I=2 pour le soluté XSEG = zeros(4,1); % Les Xj de la publi Rt = zeros(1,2); % Somme des segments du solvant (1), du soluté (2) %Xs=zeros(N,1); %Xpred=zeros(N,1); for i=1:4 for j=1:4 G(i,j) = exp( - ALPHA(i,j) * TAU(i,j)); % G de la publi end end A = load('input2.txt'); % type fichier : X Y- Y+ Z (du solvant) Xexp T Tfus deltaHfus N = size(A,1); % Nombre de points de référence sol=zeros(N,1); % initialisation de la matrice des solubilité (pour allocation) solub=zeros(N,1); % initialisation de la matrice des différences de logarithme (pour allocation) % for p=1:N, % Grosse boucle sur plusieurs solvants % Rt(1) = A(p,1)+A(p,2)+A(p,3)+A(p,4); % Le "1" d�signe le soluté. Ici, on calcule la somme des segment du solvant Rt(2) = X(1)+X(2)+X(3)+X(4); % Le "2" désigne le solvant. L� la somme des segments du soluté Xs = A(p,5); % on affecte la valeur initiale de la solubilité (expérimentale de préférence) Xpred = 0; % initialisation de la solubilité prédite err=abs(Xs-Xpred); iter=-1; while err > 10^-8, % boucle sur erreur entre les solubilit� pr�dite par le mod�le et l'�quation d'�quilibre iter=iter+1; for j=1:4 % calcul des valeurs de XSPEC et XSEG XSPEC(j,1) = 0; XSPEC(j,2) = 0; XSEG(j) = 0; XSPEC(j,1) = A(p,j)/Rt(1); XSPEC(j,2) = X(j)/Rt(2); XSEG(j) = ((A(p,j)*(1-Xs)+X(j)*Xs))/(Rt(1)*(1-Xs)+Rt(2)*Xs); end lngamma=lngammac(Xs,Rt)+lngammar(XSEG,XSPEC,G,TAU,X); Xpred = recalcul(A(p,6),A(p,7),A(p,8),lngamma); %err=abs(Xs-Xpred); %Xs=Xpred; Xpred = (2*Xs+Xpred)/3; err=abs((Xpred-Xs)/(Xs+1e-16)); Xs=Xpred; end %iter; % nombre d'iterations pour arriver au resultat disp(lngammac(Xs,Rt)); %affichage contribution entropique disp(lngammar(XSEG,XSPEC,G,TAU,X)); %affichage contribution enthalpique sol(p)=Xs; %solub(p)=(log(A(p,5))-log(sol(p)))/sqrt(N); end % Cas où on fait deja la somme des différences %rms=0; %for h=1:N % rms=rms+solub(h); %end %sol; %y=sqrt((rms)/N); % sortie somme racine carrée %y= solub; % sortie lsqnonlin y = sol; %sortie solubilité function [y]=lngammac(Xs,Rt) %calcul de la partie combinatoire PHI(1)= Rt(1)*(1-Xs)/(Rt(1)*(1-Xs)+Rt(2)*Xs); PHI(2)= Rt(2)*Xs/(Rt(1)*(1-Xs)+Rt(2)*Xs); SOM = PHI(1)/Rt(1) + PHI(2)/Rt(2); y = log( PHI(2) / Xs ) + 1 - Rt(2) * SOM; function [y] = lngammar(XSEG,XSPEC,G,TAU,X) %calcul de la partie r�siduelle dummy=0; for m=1:4 LNSACTERM1TOP = 0; LNSACTERM1BOT = 0; LNSACTERM2 = 0; LNSACITERM1TOP = 0; LNSACITERM1BOT = 0; LNSACITERM2 = 0; for j=1:4 LNSACTERM2BOT = 0; LNSACTERM4TOP = 0; LNSACTERM4BOT = 0; LNSACITERM2BOT = 0; LNSACITERM4TOP = 0; LNSACITERM4BOT = 0; for k=1:4 LNSACTERM2BOT = LNSACTERM2BOT + XSEG(k)*G(k,j); LNSACTERM4TOP = LNSACTERM4TOP + XSEG(k)*G(k,j)*TAU(k,j); LNSACTERM4BOT = LNSACTERM4BOT + XSEG(k)*G(k,j); LNSACITERM2BOT = LNSACITERM2BOT + XSPEC(k,2)*G(k,j); LNSACITERM4TOP = LNSACITERM4TOP + XSPEC(k,2)*G(k,j)*TAU(k,j); LNSACITERM4BOT = LNSACITERM4BOT + XSPEC(k,2)*G(k,j); end LNSACTERM1TOP = LNSACTERM1TOP + XSEG(j)*G(j,m)*TAU(j,m); LNSACTERM1BOT = LNSACTERM1BOT + XSEG(j)*G(j,m); LNSACTERM2 = LNSACTERM2 + XSEG(j) * G(m,j)/LNSACTERM2BOT * ( TAU(m,j) - LNSACTERM4TOP / LNSACTERM4BOT ); LNSACITERM1TOP = LNSACITERM1TOP + XSPEC(j,2)*G(j,m)*TAU(j,m); LNSACITERM1BOT = LNSACITERM1BOT + XSPEC(j,2)*G(j,m); LNSACITERM2 = LNSACITERM2 + XSPEC(j,2)*G(m,j)/LNSACITERM2BOT * ( TAU(m,j) - LNSACITERM4TOP / LNSACITERM4BOT); end LNSAC = LNSACTERM1TOP / LNSACTERM1BOT + LNSACTERM2; LNSACI = LNSACITERM1TOP / LNSACITERM1BOT + LNSACITERM2; dummy = dummy + X(m) * ( LNSAC - LNSACI ); y = dummy; end function [y] = recalcul(T,Tfus,Hfus,lngamma) y=exp(Hfus/8.314*(1/Tfus-1/T)-lngamma); % recalcul de la solubilité avec l'équation d'équilibre II.2/NRTL-SAC ternaire FORTRAN PROGRAM NRTL_SAC !****************************************************************************** ! GENERAL COMMENTS: ! Author: Adel Ghaderi ! Supervisor: Eric Mullins ! Advisor: Y. A. Liu ! Institution: Virginia Polytechnic Institute and State University ! Date: October 27,2006 ! Code is setup to handle pure solvent systems for solubility of single solute ! for subscripts referring to the segments ! 1 refers to X ! 2 refers to Y- ! 3 refers to Y+ ! 4 refers to Z ! Program is coded for 4 segments; simply changing value of NUMSEG is not ! sufficient to vary number of segments. !****************************************************************************** ! REFERENCES: ! Chen, C.C.; Song, Y.; "Solubility modeling with a nonrandom two-liquid segment activity coefficient ! model," IECR (2004), v.43, no.26, 8354-8362. !****************************************************************************** IMPLICIT NONE INTEGER, PARAMETER :: SINGLE = 4 INTEGER, PARAMETER :: DOUBLE = 8 INTEGER :: i, j, k, m, n, y, z, ITER, ITERcount, NUMSPEC, NUMTEMP, NUMSYS, istat INTEGER, PARAMETER :: NUMSEG = 4 !number of segments; 4 (X, Y-, Y+, Z) as DEFINED by C. C. Chen INTEGER, PARAMETER :: CHARLEN = 100 !maximum allowable length name of the system REAL(KIND=DOUBLE) :: LNSACTERM1TOP, LNSACTERM1BOT, LNSACTERM2 REAL(KIND=DOUBLE) :: LNSACITERM1TOP, LNSACITERM1BOT, LNSACITERM2 REAL(KIND=DOUBLE) :: LNSACTERM2BOT, LNSACTERM4TOP, LNSACTERM4BOT REAL(KIND=DOUBLE) :: LNSACITERM2BOT, LNSACITERM4TOP, LNSACITERM4BOT REAL(KIND=DOUBLE) :: LNSAC, LNSACI, dummy, XOLD, MFSUM, RELTOL, CHANGE, LNRE, SOM, Cp REAL(KIND=DOUBLE), ALLOCATABLE, DIMENSION(:) :: T, XSEG, LNGAMMAR, LNGAMMAC, PHI, dHfus, TM, Rt, X, XORIG !, Cp REAL(KIND=DOUBLE), ALLOCATABLE, DIMENSION(:,:) :: TAU, ALPHA, G, R, XSPEC, GAMMA CHARACTER(CHARLEN), ALLOCATABLE, DIMENSION(:) :: TITLE CHARACTER(CHARLEN), ALLOCATABLE, DIMENSION(:,:) :: NAME !------------------------------------------------------------------------------------ ! Initialize Tau and Alpha values, from C.C. Chen, et al. !------------------------------------------------------------------------------------ ALLOCATE( TAU(NUMSEG,NUMSEG), ALPHA(NUMSEG,NUMSEG), G(NUMSEG,NUMSEG) ) DO i = 1, 4 ! 4 segment species (X, Y-, Y+, Z) DO j = 1, 4 TAU(i,j) = 0 ALPHA(i,j) = 0 ! note ALPHA(i,j) = ALPHA(j,i) END DO END DO TAU(1,2) = 1.643 ! 1 denotes X (hydrophobic segment) TAU(2,1) = 1.834 ! 2 denotes Y- (polar attractive segment) ALPHA(1,2) =0.200; ALPHA(2,1) = .200 TAU(1,4) = 6.547 ! 4 denotes Z (hydrophilic segment) TAU(4,1) = 10.949 ALPHA(1,4) =0.200; ALPHA(4,1) = .200 TAU(2,4) = -2.000 TAU(4,2) = 1.787 ALPHA(2,4) =0.300; ALPHA(4,2) = .300 TAU(3,4) = 2.000 ! 3 denotes Y+ (polar repulsive segment) TAU(4,3) = 1.787 ALPHA(3,4) =0.300; ALPHA(4,3) = .300 TAU(1,3) = 1.643 TAU(3,1) = 1.834 ALPHA(1,3) =0.200; ALPHA(3,1) = .200 !------------------------------------------------------------------------------------ ! Calculate values for G DO i = 1, 4 DO j = 1, 4 G(i,j) = exp( - ALPHA(i,j) * TAU(i,j) ) END DO END DO ! Generate input file OPEN(UNIT = 13, FILE = "base_input_ternaire_cp.txt", & STATUS="OLD", ACTION="READ", POSITION="REWIND") OPEN(UNIT=14,FILE="input_intermediaire_ternaire_cp.txt",STATUS="REPLACE", ACTION="WRITE") ! The input data of interest is for pure component solubility NUMSPEC = 3 ! For simplicity, every point is treated as a single system (this is usually the case anyway) NUMTEMP = 1 119 FORMAT( I4 ) 120 FORMAT( I1, 5X, I1 ) 121 FORMAT( A15 ) 122 FORMAT( F10.2,2X,F10.2,2X,F7.2,2X,F12.10,2X,F12.10,2X,F12.10,2X,F12.10 ) 123 FORMAT( F8.2, 2X, F17.15 ) 124 FORMAT( F10.2,2X,I1,2X,I1,2X,F12.10,2X,F12.10,2X,F12.10,2X,F12.10 ) READ(13,119) NUMSYS !NUMSYS = 1 WRITE(14,119) NUMSYS ALLOCATE( R(12,NUMSYS), dHfus(NUMSYS), TM(NUMSYS), T(NUMSYS), X(NUMSPEC), NAME(NUMSYS,NUMSPEC), STAT = istat ) WRITE(*,*) istat , " ALLOCATION SUCCESSFUL IF 0" DO i = 1, NUMSYS ! INPUT FILE FORMAT: SOLUTE NAME, HEAT OF FUSION, MELTING TEMPERATURE, SEGMENT PARAMETERS (X,Y-, ! Y+,Z),SOLVENT NAME, SEGMENT PARAMETERS (X,Y-,Y+,Z), SYSTEM TEMPERATURE, EXPERIMENTAL SOLUBILITY READ(13,*) NAME(i,1), dHfus(i), TM(i), Cp, R(1,i), R(2,i), R(3,i), R(4,i), NAME(i,2), R(5,i), R(6,i), & R(7,i), R(8,i), NAME(i,3), R(9,i), R(10,i), R(11,i), R(12,i), T(i), X(2) X(3) = 0.0000001D0 WRITE(14,120) NUMSPEC, NUMTEMP WRITE(14,121) NAME(i,2) WRITE(14,121) NAME(i,3) WRITE(14,121) NAME(i,1) WRITE(14,124) Cp, 0, 0, R(5,i), R(6,i), R(7,i), R(8,i) WRITE(14,124) Cp, 0, 0, R(9,i), R(10,i), R(11,i), R(12,i) WRITE(14,122) Cp, dHfus(i), TM(i), R(1,i), R(2,i), R(3,i), R(4,i) WRITE(14,123) T(i), X(2) END DO DEALLOCATE( R, dHfus, TM, T, X, NAME, STAT = istat ) WRITE(*,*) istat, " DEALLOCATION SUCCESSFUL IF 0" CLOSE(13) CLOSE(14) !------------------------------------------------------------------------------------ ! Open input and output files OPEN(UNIT = 11, FILE = "input_intermediaire_ternaire_cp.txt",STATUS="OLD", ACTION="READ", POSITION="REWIND") READ(11,*) NUMSYS OPEN(UNIT=12, FILE = "output_ternaire_cp.txt",& STATUS="REPLACE", ACTION="WRITE") OPEN(UNIT=15, FILE = "resultat_fraction_cp.txt", STATUS = "REPLACE", ACTION = "WRITE") WRITE(12,115) "ITER", "SOLUTE", "SOLVENT", "TEMP(K)", "GAMMA-3", "X3", "Rel.-%E-(X2)" WRITE(12,115) "====", "======", "=======", "=======", "=======", "==", "============" 115 FORMAT(2X, A4, 3X, A6, 6X, A7, 4X, A7, 4X, A7, 8X, A2, 8X, A12) ! Iterate over all different multicomponent systems !------------------------------------------------------------------------------------ ALLOCATE( NAME(NUMSPEC,1), TITLE(NUMSYS), dHfus(NUMSPEC), TM(NUMSPEC), & T(1), XSEG(NUMSEG), LNGAMMAR(NUMSPEC), LNGAMMAC(NUMSPEC), PHI(NUMSPEC), & R(NUMSEG,NUMSPEC), XSPEC(NUMSEG,NUMSPEC), XORIG(NUMSPEC), X(NUMSPEC), & GAMMA(NUMSEG,NUMSPEC), Rt(NUMSPEC), STAT = istat ) ! WRITE(*,*) istat , " ALLOCATION SUCCESSFUL IF 0" DO z = 1, NUMSYS !number of system components and number of data points; used in memory allocation READ(11,*) NUMSPEC, NUMTEMP WRITE(*,*) NUMSPEC, NUMTEMP ! Collect VLE data from file DO i = 1, NUMSPEC READ(11,*) NAME(i,1) END DO DO i = 1, NUMSPEC READ(11,*) Cp , dHfus(i), TM(i), R(1,i), R(2,i), R(3,i), R(4,i) ! values for X, Y-, Y+, and Z for species i END DO READ(11,*) T(1), XORIG(2) ! X(1) = XORIG(1) ! Pour utiliser deux compositions dans le fichier d'entrée, utiliser le progamme "deux_compo..." X(1) = 1- XORIG(2) X(2) = XORIG(2) X(3) = 0.00000001D0 ! Generate XSPEC, denoted Xj,I in C. C. Chen (note: XSPEC is independent of system composition) et somme R(i,j) !------------------------------------------------------------------------------------ DO i = 1, NUMSPEC dummy = 0 DO j = 1, 4 XSPEC(j,i) = 0 Rt(i)=0 DO k = 1, 4 Rt(i) = Rt(i) + R(k,i) END DO XSPEC(j,i) = R(j,i) / Rt(i) END DO END DO ! Boucle sur la température DO n = 1, NUMTEMP ! Début de la boucle de détarmination de la solubilité DO ITER = 1, 200 !if (z==2) then !read(*,*) j !endif ! Generate XSEG, denoted Xj in C. C. Chen !------------------------------------------------------------------------------------ DO i = 1, 4 dummy = 0 XSEG(i) = 0 DO j = 1, NUMSPEC dummy = dummy + X(j)*Rt(j) XSEG(i) = XSEG(i) + X(j)*R(i,j) END DO XSEG(i) = XSEG(i) / dummy END DO ! calculate over all species DO i = 1, NUMSPEC ! Calculate residual term, LNGAMMAR, from NRTL-SAC model !------------------------------------------------------------------------------------ LNGAMMAR(i) = 0 DO m = 1, 4 LNSACTERM1TOP = 0 LNSACTERM1BOT = 0 LNSACTERM2 = 0 LNSACITERM1TOP = 0 LNSACITERM1BOT = 0 LNSACITERM2 = 0 DO j = 1, 4 !j is same as j, k, and m' for terms 1 and 2 LNSACTERM2BOT = 0 LNSACTERM4TOP = 0 LNSACTERM4BOT = 0 LNSACITERM2BOT = 0 LNSACITERM4TOP = 0 LNSACITERM4BOT = 0 DO k = 1, 4 !k is same as k, j, and k for terms 2 and 4 LNSACTERM2BOT = LNSACTERM2BOT + XSEG(k)*G(k,j) LNSACTERM4TOP = LNSACTERM4TOP + XSEG(k)*G(k,j)*TAU(k,j) LNSACTERM4BOT = LNSACTERM4BOT + XSEG(k)*G(k,j) ! note same as LNSACTERM2BOT LNSACITERM2BOT = LNSACITERM2BOT + XSPEC(k,i)*G(k,j) LNSACITERM4TOP = LNSACITERM4TOP + XSPEC(k,i)*G(k,j)*TAU(k,j) LNSACITERM4BOT = LNSACITERM4BOT + XSPEC(k,i)*G(k,j) ! note same as LNSACITERM2BOT END DO ! end iteration of k !note that G(m,j) in LNSACTERM2 and LNSACITERM2 have indicies reversed !from those of other G values, but this is confirmed by lit (C.C. Chen) LNSACTERM1TOP = LNSACTERM1TOP + XSEG(j)*G(j,m)*TAU(j,m) LNSACTERM1BOT = LNSACTERM1BOT + XSEG(j)*G(j,m) LNSACTERM2 = LNSACTERM2 + XSEG(j) * G(m,j) / & LNSACTERM2BOT * ( TAU(m,j) - LNSACTERM4TOP / LNSACTERM4BOT ) LNSACITERM1TOP = LNSACITERM1TOP + XSPEC(j,i)*G(j,m)*TAU(j,m) LNSACITERM1BOT = LNSACITERM1BOT + XSPEC(j,i)*G(j,m) LNSACITERM2 = LNSACITERM2 + XSPEC(j,i) * G(m,j) / & LNSACITERM2BOT * ( TAU(m,j) - LNSACITERM4TOP / LNSACITERM4BOT ) END DO ! end iteration of j LNSAC = LNSACTERM1TOP / LNSACTERM1BOT + LNSACTERM2 LNSACI = LNSACITERM1TOP / LNSACITERM1BOT + LNSACITERM2 !if (i==2) then ! write(*,*) 'm=',m,'lnac =', LNSAC, 'lnaci =', LNSACI !end if LNGAMMAR(i) = LNGAMMAR(i) + R(m,i) * ( LNSAC - LNSACI ) END DO ! end iteration of m !------------------------------------------------------------------------------------ ! Calculate combinatorial term, LNGAMMAC, from Flory-Huggins, as done by C. C. Chen !------------------------------------------------------------------------------------ ! calcul de PHI(k) DO k = 1, NUMSPEC PHI(k) = 0 DO j = 1, NUMSPEC PHI(k) = PHI(k) + Rt(j)*X(j) END DO PHI(k) = Rt(k)*X(k)/PHI(K) END DO ! calcul de SOM = somme sur k de [PHI(k)/Rt(k)] SOM = 0 DO j = 1, NUMSPEC SOM = SOM + PHI(j)/Rt(j) END DO ! calcul de la partie combinatoire, model as shown by C. C. Chen LNGAMMAC(i) = LOG( PHI(i) / X(i) ) + 1 - Rt(i) * SOM !write(*,*) 'SOM =', SOM, 'ln gamma résid', i, '=', LNGAMMAC(i) !------------------------------------------------------------------------------------ ! Calculate overall activity coefficient, GAMMA, and predicted Y, YPRED !------------------------------------------------------------------------------------ GAMMA(n,i) = EXP( LNGAMMAR(i) + LNGAMMAC(i) ) !------------------------------------------------------------------------------------ END DO ! end iteration of i ! XOLD = X(NUMSPEC) ! RECALCULATE SOLUTE MOLE FRACTION FROM DERIVED EQUATION X(NUMSPEC)=EXP(dHfus(NUMSPEC)/(TM(NUMSPEC)*8.314)*(1-TM(NUMSPEC)/T(1))- & Cp/8.314*(LOG(TM(NUMSPEC)/T(n))-TM(NUMSPEC)/T(n)+1) -LOG(GAMMA(n,NUMSPEC))) X(NUMSPEC) = ( 2.0D0*XOLD + X(NUMSPEC) ) / 3.0D0 !CHECK FOR NEGATIVE MOLE FRACTIONS IF (X(NUMSPEC) >= 1.0) THEN WRITE (*,*) "SOLUTE MOLEFRACTION CALCULATED AS GREATER THAN 1.0" WRITE (*,*) X(NUMSPEC) END IF !NORMALIZE MOLE FRACTIONS KEEPING SOLVENT RATIO (X1:X2) CONSTANT MFSUM = 0 DO y = 1, NUMSPEC MFSUM = MFSUM + X(y) END DO DO y = 1, NUMSPEC X(y) = X(y) / MFSUM END DO ! CALCULATE THE RELATIVE CHANGE OF THE SOLUTE MOLE FRACTION CHANGE = ABS( ( X(NUMSPEC) - XOLD)/( XOLD + 1.0D-16 ) ) !WRITE (*,*) " RELATIVE CHANGE:", CHANGE !WRITE(*,*) ITER, CHANGE, X(1),X(2),X(3) ! COMPARE RELATIVE TOLERANCE RELTOL = 1.E-8 IF (CHANGE < RELTOL) THEN ITERcount = ITER GOTO 200 END IF !END ITERATIVE LOOP END DO ! END OF ITERATION UPON ITER 200 dummy = ABS( X(NUMSPEC) - XORIG(NUMSPEC) ) / XORIG(NUMSPEC) * 100 ! relative error ! Ecriture des résultats dans le fichier 12 et dans le shell WRITE(12,111) ITER, NAME(3,1), NAME(1,1), NAME(2,1), T(1), GAMMA(n,NUMSPEC), XORIG(NUMSPEC),& X(NUMSPEC) ! Pure Solvent WRITE(*,111) ITER, NAME(3,1), NAME(1,1), NAME(2,1), T(1), GAMMA(n,NUMSPEC),& XORIG(NUMSPEC), X(NUMSPEC) 111 FORMAT(2X,I3,4X,A9,1X,A10,1X,A10,1X,F8.2,2X,F12.8,2X,F8.5,5X,F8.5) dummy = 0 WRITE(*,*) "solubility =", X(NUMSPEC) WRITE(15,*) X(NUMSPEC) END DO ! end iteration of n END DO ! end iteration of z DEALLOCATE( NAME, TITLE, dHfus, TM, T, XSEG, LNGAMMAR, LNGAMMAC, PHI, R, XSPEC, XORIG, X,& GAMMA, Rt, TAU, ALPHA, G, STAT = istat) CLOSE(11) CLOSE(12) CLOSE(15) WRITE(*,*) "DONE!" END PROGRAM NRTL_SAC