c Compile: gfortran NVT_paper.f c Run: ./a.out c Inputs: 3 (for sphere), 2 (for cylinder), 1 (for planar slab) c Outputs: c - Number of crystalline particles (for R=Rs), c - Change in F between homogeneous liquid and heterogeneous system c - Pressure of the solid c - Density of the solid c - Pressure of the liquid c - Density of the liquid c - Radius at the surface of tension, Rs c - As*gamma_s/3 (the Gibbs free energy barrier to nucleation when F is a minimum c - Inverse of the radius (for R=Rs) ******************************************************************************************** implicit real*8 (a-h,o-z) pi = dacos(-1.00d0) c ****************************************************************************************** c The program will use reduced units of HS c p=p/(kT/sigma**3) rho=(N/V)*sigma**3 c Sigma is used as unit of length: c Length in sigma c Area in sigma**2 c Volume in sigma**3 c Interfacial free energy in kT/sigma**2 c Free energies (F and G) and chemical potetials in kT units c where F is Helmholtz free energy and G is Gibbs free energy open(unit=22,file="Nsol_DF_ps_rhos_pl_rhol_Gc.dat", 1 status="unknown") c write(22,*) c 1 "Nsol_Rs", " ", c 1 "Fhet-Fhom", " ", c 1 "psol", " ", c 1 "rhos", " ", c 1 "pliq", " ", c 1 "rhol", " ", c 1 "Rs", " ", c 1 "As*gamma_s/3", " ", c 1 "1/Rs" c ****************************************************************************************** c The number of particles (N) and volume (V) of the ten systems considered in this work, c previously studied in P.Montero et al., J.Phys.Chem.C, (2020) are: c System I V=10686.4 xN=10540. c System II c V=20195.5 c xN=19779 c System III c V=20195.5 c xN=19829 c System IV c V=49599.9 c xN=48207. c System V c V=49599.9 c xN=48357 c System VI c V=108265.2 c xN=104675 c System VII c V=66900.1 c xN=65383 c System VIII c V=108265.2 c xN=105475 c System IX c V=108265.2 c xN=105875 c System X c V=887000.0 c xN=853712 c Side of the simulation box xbox = V**(1./3.) c ****************************************************************************************** c Defining the geometry of the solid-liquid interface write(6,*) " Input geometry of the solid cluster " write(6,*) " 1 Plane 2 Cylinder 3 Sphere " read(5,*) igeometry c ****************************************************************************************** c Equation of state of the fluid phase P = 141.04*rho**2 - 212.12*rho + 86.29 al = 141.04 bl = -212.12 cl = 86.29 c Equation of state of the solid phase P = 162.96*rho**2 - 299.42*rho + 146.8 as = 162.96 bs = -299.42 cs = 146.8 c ****************************************************************************************** c Interfacial free energy of the planar interface at coexistence c (average value of several planes) g0 = 0.576*1.02 c Tolman length tolman = -0.41 c ****************************************************************************************** c Parameters that can be tunned for accuracy c Width of the grid for radius deltaR = 0.005 c Width of the grid for pressure deltaP = 0.00125 c Number of pressure points used in solving the equation N = N_solid + N_liquid npmax = 10./deltaP c Initial pressure p0 = 8.5 c Maximum number of points used for the radius nrmax = xbox/2./deltaR c ****************************************************************************************** c Thermodynamic properties at coexistence (without curvature, hereafter simply coexistence) c Coexistence pressure p_cx = 11.648 !espinosa et al. xmu_cx = 0.0 ! coex.chem.pot. set to zero arbitrarily clc = cl-p_cx c Sinde p is a quadratic function of rho, to determine c the value of rho for a certain p we solve the quadratic equation rhol_cx = (-bl+dsqrt(bl**2-4.*al*clc))/(2.*al) ! Dens.liq.coex. Fl_cx = xmu_cx - p_cx/rhol_cx ! F of the liquid (per particle) csc = cs-p_cx rhos_cx = (-bs+dsqrt(bs**2-4.*as*csc))/(2.*as) ! Den.sol.coex. Fs_cx = xmu_cx - p_cx/rhos_cx ! F of the solid (per particle) c ****************************************************************************************** c Properties of the homogeneous fluid in the simulation box having N and V conditions c Density rhol_bk = xN/V c Packing fraction phi_bk = (pi/6.)*rhol_bk c Pressure pl_bk = al*rhol_bk**2 + bl*rhol_bk + cl c Evaluating the change from coexistence in F for the liquid using thermodynamic integration c (see Eq (5) from Vega et al. J. Phys. Condens. Matter (2008)) rho1 = rhol_cx rho2 = rhol_bk delta = al*(rho2-rho1)+bl*dlog(rho2/rho1)+cl*(1./rho1-1./rho2) c where delta is the change in F for a liquid that goes from coexistence to pl_bk c Helmholtz free energy of the homogeneous liquid Fl_bk = Fl_cx + delta c Chemical potential of the homogeneous liquid xmu_bk = Fl_bk + pl_bk/rhol_bk c ****************************************************************************************** c Some geometry-dependent parameters: xg = igeometry if (igeometry.eq.3) then alpha = (36.*pi)**(1.d0/3.d0) expo = 2.d0/3.d0 nmax = (xbox/2.)/deltaR elseif (igeometry.eq.2) then alpha = 2.*dsqrt(pi*xbox) expo = 1.d0/2.d0 nmax = (xbox/2.)/deltaR elseif (igeometry.eq.1) then alpha = 2.*(xbox)**2 expo = 0. nmax = xbox/deltaR endif c ****************************************************************************************** c Loop over different values of the radius c Rs is the radius at the surface of tension c Re is the radius at the equimolar dividing surface do i = 750,nrmax Rs = i*deltaR Re = Rs + tolman + 2.613/Rs c Interfacial free energy at the surface of tension via Tolman's equation g_s = g0*(1.-(xg-1.)*tolman/Rs) c Interfacial free energy at the equimolar dividing surface g_e = g_s*((Rs**xg + (xg-1.)*Re**xg)/xg/Rs/Re**(xg-1.)) c Flag variables xdif_old = 0. iget = 0 c Loop over possible trial values of the solid pressure for a given Rs do j = 1,npmax psol = p0 + (j-1)*deltaP 333 if (iget.eq.1) then !N = N_solid + N_liquid was solved, thus: psol = psol_final endif c Pressure of the fluid phase according to Young-Laplace equation pliq = psol-(xg-1.)*g_s/Rs c Properties of the liquid at pliq cll = cl-pliq xl = bl**2-4.*al*cll if (xl.gt.0.) then rhol = (-bl+dsqrt(xl))/(2.*al) !Density of the liquid rho1 = rhol_cx rho2 = rhol delta = al*(rho2-rho1)+bl*dlog(rho2/rho1)+ 1 cl*(1./rho1-1./rho2) c where delta is now the change in F for a liquid that goes from coexistence to pliq Fliq = Fl_cx+delta !F of the liquid xmu_liq = Fliq+pliq/rhol !chemical potential of the liquid c Properties of the solid phase at psol css = cs-psol xs = bs**2-4.*as*css if (xs.gt.0.) then rhos = (-bs+dsqrt(xs))/(2.*as) !Density of solid rho1 = rhos_cx !Density of the solid at coexistence rho2 = rhos delta = as*(rho2-rho1)+bs*dlog(rho2/rho1)+ 1 cs*(1./rho1-1./rho2) c where delta is now the change in F for a solid that goes from coexistence to psol Fsol = Fs_cx + delta !F of the solid xmu_sol = Fsol + psol/rhos ! Chemical potential of the solid c Determining volumes and number of particles according to the geometry c Volumes and number of particles using both Re and Rs: !Solid if (igeometry.eq.3) then Vs_e = (4./3.)*pi*Re**3 Vs_s = (4./3.)*pi*Rs**3 elseif (igeometry.eq.2) then Vs_e = (pi*Re**2)*xbox Vs_s = (pi*Rs**2)*xbox elseif (igeometry.eq.1) then Vs_e = Re*xbox**2 Vs_s = Rs*xbox**2 endif !Liquid Vl_e = V-Vs_e Vl_s = V-Vs_s xns_e = Vs_e*rhos xns_s = Vs_s*rhos xnl_e = Vl_e*rhol xnl_s = Vl_s*rhol c Evaluating the number of particles in the system using the trial value of Re xtry = Vs_e*rhos + Vl_e*rhol c Difference between actual number of particles and the calculated one xdif = xN-xtry c At the true value of psol this difference should be zero c When the difference between actual and calculated number of particles change its c sign the value at which is zero should be between them c iget is a flag for convergence (0 no previous convergence, 1 yes previous convergence) xtest = xdif_old*xdif if (xtest.lt.0. .and. iget .eq. 0) then c Linear interpolation to determine the pressure at which the c difference between actual and calculated number of particles is zero c (Nexacta - Ncalculada) = 0 c slope = (y''-y')/(x''-x'); y2 = 0 -> x2 = x1 - y1/slope slope = (xdif-xdif_old)/deltaP psol_final=(psol-deltap)-xdif_old/slope iget = 1 !Convergence obtained go to 333!Go to Final iteration endif if (iget .eq. 1) then ! Determining final properties c Area of the interface at Re As_e = alpha*Vs_e**expo c Helhmoltz free energy F per particle of the heterogeneous system F_het = (xnl_e*Fliq+xns_e*Fsol+As_e*g_e)/xN c Helhmoltz free energy F per particle of the homogeneous system F_hom = Fl_bk dF = xN*(F_het - F_hom) !Difference in F c Area of the interface at Rs As_s = alpha*Vs_s**expo c At the free energy minimum in F, 1/3 Area * gamma_s gives Delta G_c c the free energy for the formation of the critical cluster at c constant pressure, useful for nucleation studies deltagc=1./3.*As_s*g_s c Final output Nsol(Rs),Delta F, psol, rhosol, pliq, rholiq, Rs, c 1/3 interfacial free energy which is Delta G_c at the minimum in F write(22,741) Vs_s*rhos, 1 dF, 1 psol, 1 rhos, 1 pliq, 1 rhol, 1 Rs, 1 deltagc, 1 1./Rs go to 20 ! psol for Rs found. Go to compute the next Rs endif endif xdif_old=xdif endif end do ! Start another trial psol 20 end do write(6,*) "END of program" 741 format(9(f16.6,2x)) end