Calculating surface tension¶

In this tutorial we demonstrate how to calculate the surface tension between water and hexadecane.

The first place to start is the definition of the surface tension $γ$

\frac{d F}{d A} \equiv 2 γ

where $F$ is the Helmholtz free energy and $A$ is the area of the interface. Our approach will be to calculate the $d F / d A$ derivative using a finite difference.

The Helmholtz free energy can be conveniently found using SCFT by locating the saddle point fields $w^{*}$ :

β F = - \ln Z_{c} \approx - \ln (e^{- H [w^{*}]}) = H [w^{*}] .

To begin, lets consider the units in OpenFTS. OpenFTS provides two options to specify units.

Rg: units are given in terms of a reference radius of gyration $R_{g, r e f} = b_{r e f} \sqrt{N_{r e f} / 6}$ (continuous chains) or $R_{g, r e f} = b_{r e f} \sqrt{(N_{r e f} - 1) / 6}$ (discrete chains)
b : units are given in terms of a reference statistical segement length $b_{r e f}$

We will use Rg units and will use the convention that $R_{g, r e f} = 1$ nanometer.

[1]:

units = 'Rg'

In our model, we will represent hexadecane by a 16 bead polymer where each bead represents a CH2 monomer (we ignore end effects). We will represent water by a single bead.

[2]:

N = 16

The statistical segment length of hexadecane is $b_{C} = 0.33565$ nm. Since water is a single bead, the statistical segment length will not matter so we will simply set $b_{W} = 1$ nm.

[3]:

b_C = 0.137 # units of bref = 2.45 nm
b_W = 1

For simplicity, we will also set the smearing length for hexadecane and water to $a_{C} = a_{W} = 1$ nm.

[4]:

a_C = a_W = 1 # nm

Finally, we need to specify the interaction parameters between the C and W monomers.

[5]:

chi = 2.460621 # this was found using COSMO-RS
zeta = 100 # Arbitrarily specified

Now that all parameters have been specified, we now define a function that runs SCFT for a specified $L_{x}$ (x cell length), $A_{y z}$ (yz cross sectional area) and $d A$ (perturbation in cross sectional area):

[ ]:

import openfts
import numpy as np

def run(Lx, Ayz, dA=0):

    # perturb yz area while keeping total volume constant
    Ly = Lz = np.sqrt(Ayz + dA)
    Lx = Lx*Ayz/(Ly*Lz)

    fts = openfts.OpenFTS()

    fts.cell = openfts.Cell(cell_scale=1.0, cell_lengths=[Lx,Ly,Lz], dimension=3,
                            tilt_factors=[0, 0, 0], length_unit=units)
    fts.field_layout = openfts.FieldLayout(npw=[100, 20, 20], random_seed=777)

    scft = openfts.driver.SCFT(dt=1.0, nsteps=100000, output_freq=1000, stop_tolerance=1e-5)
    scft.field_updater = openfts.field_updater.EMPEC(update_order='simultaneous',
                                                     adaptive_timestep=False, lambdas=[0.02, 0.1])
    fts.driver = scft

    chiMultiSpecies = openfts.model.ChiMultiSpecies(Nref=2, bref=1.0, chi_array=[chi], zeta=zeta, C=5)
    chiMultiSpecies.init_fields = openfts.init_field.FromFile('fields.in', field_type='exchange')
    fts.model = chiMultiSpecies

    fts.molecules.append(openfts.molecule.PolymerLinearDiscrete(monomer_sequence='C'*N, nbeads=N, volume_frac=0.5))
    fts.molecules.append(openfts.molecule.Bead(bead_species='W', volume_frac=0.5))
    fts.operators.append(openfts.operator.Hamiltonian(units='beta_Rg3_over_V'))
    fts.output = openfts.Output()
    fts.species.append(openfts.Species(label='C', stat_segment_length=b_C,
                       smearing_length=a_C, smearing_length_units=units))
    fts.species.append(openfts.Species(label='W', stat_segment_length=b_W,
                       smearing_length=a_W, smearing_length_units=units))
    fts.run()

    H = openfts.get_scalar_operator('operators.dat', 'H.real', 'final')
    return H

Now lets run SCFT for a specified Lx, Ayz and dA=0 and save the result in H1

[7]:

Lx = 100 #nm
Ayz = 400 #nm
H1 = run(Lx, Ayz, dA=0.0)
print(H1)

================================================================
                           OpenFTS

OpenFTS version: heads/develop f48f964 DIRTY
FieldLib version: heads/master c3fc5fd CLEAN
FieldLib precision: double
Licensee: Joshua Lequieu | jl4347@drexel.edu | Expiration: 12/31/25
================================================================

Setup Species
  * label = C
  * stat_segment_length (b) = 0.137
  * smearing_length (a) = 1
  * smearing_length_units = Rg
  * charge = 0.000000
Setup Species
  * label = W
  * stat_segment_length (b) = 1
  * smearing_length (a) = 1
  * smearing_length_units = Rg
  * charge = 0.000000
Initialize Cell
  * dim = 3
  * length unit = Rg
  * cell_lengths = [100.000000 20.000000 20.000000 ]
  * tilt_factors = [0.000000 0.000000 0.000000 ]
  * cell_tensor:
      [ 100,   0,   0 ]
      [   0,  20,   0 ]
      [   0,   0,  20 ]
Initialize FieldLayout
  * npw = [100 20 20 ]
  * random_seed = 777
  * backend = NOT_SPECIFIED
  * backend = CUDA (chosen by fieldlib)
  * Execution device: [GPU 0] NVIDIA GeForce RTX 2080 Ti (ComputeCapability 7.5)
Setup MoleculePolymerLinear
  * chain_type = discrete
  * label = mol0
Setup MoleculeBead
  * label = mol1
  * volume_frac = 0
  * bead_species = W
Setup ModelChiMultiSpecies
  * bref = 1
  * Nref = 2
  * C = 5
  * rho0 = 146.969
  * compressibility_type = Helfand
  * inverse_zetaN = 0.005
  * zeta = 100
  * chiN_matrix =
[ [       0, 4.92124 ],
  [ 4.92124,       0 ]]
  * X_matrix =
[ [ 14.1421, 14.4901 ],
  [ 14.4901, 14.1421 ]]
  * O_matrix (cols are eigenvectors) =
[ [ -0.707107,  0.707107 ],
  [  0.707107,  0.707107 ]]
  * eigenvalues = -0.347984 28.6323
  * A_matrix =
[ [ -0.707107,  0.707107 ],
  [  0.707107,  0.707107 ]]
  * nfields = 2
  * has_electrostatic_field = 0
Initialize MoleculePolymerLinear
  * Linear polymer has single chain:
    - type = discrete
    - nbeads = 16
    - discrete_type = Gaussian
    - sequence_type = monomers
    - nblocks = 1
    - block_fractions =  1.000000
    - block_species =  C
    - monomer_sequence (compact) =  C_16
    - monomer_sequence (long) = CCCCCCCCCCCCCCCC
    - species_fractions =  1 0
  * volume_frac = 0.5
  * ncopies = 12500
Initialize MoleculeBead
  * volume_frac = 0.5
  * ncopies = 200000
Initialize ModelChiMultiSpecies
  * total net_charge / (CV) = 0
  * exchange_field_names = 'mu_1' 'mu_2'
  * initfields = yes
    - reading 'exchange' fields from file 'fields.in'
Initialize DriverSCFTCL
  * is_complex_langevin = false
  * dt = 1
  * output_freq = 1000
  * output_fields_freq = 1000
  * block_size = 1000
  * field_updater = EMPEC
    - predictor_corrector = true
    - update_order = simultaneous
    - lambdas = 0.020000 0.100000
    - adaptive_timestep = false

  * field_stop_tolerance = 1e-05
Calculating Operators:
  * OperatorHamiltonian
    - type = scalar
    - averaging_type = none
    - save_history = true
    - units = beta_Rg3_over_V
Output to files:
  * scalar_operators  -> "operators.dat"
  * species_densities -> "density.dat"
                      ->  save_history = false
  * exchange_fields   -> "fields.dat"
                      ->  save_history = false
  * field_errors      -> "error.dat"

Running DriverSCFTCL
  * nsteps = 100000
  * resume = false
  * store_operator_timeseries = true

Starting Run.
Step:      1  H: -2.509370e+00 -2.609e-16i  FieldErr: 5.44e+00  TPS:  59.58
Step:   1000  H: -2.391474e+00 +1.748e-17i  FieldErr: 3.87e-04  TPS: 141.68
Step:   2000  H: -2.391465e+00 -1.947e-16i  FieldErr: 6.19e-05  TPS: 144.59
Step:   3000  H: -2.391464e+00 +5.016e-16i  FieldErr: 1.39e-05  TPS: 142.83
SCFT converged in 3230 steps to tolerance: 0.00 (fields)
* H = -2.39146 -2.20678e-16i
Run completed in 23.20 seconds (139.18 steps / sec).
-2.391464122693

Next lets run SCFT for a specified Lx, Ayz and dA=4 and save the result in H2

[8]:

dA = 4 # nm2
H2 = run(Lx, Ayz, dA=dA)
print(H2)

================================================================
                           OpenFTS

OpenFTS version: heads/develop f48f964 DIRTY
FieldLib version: heads/master c3fc5fd CLEAN
FieldLib precision: double
Licensee: Joshua Lequieu | jl4347@drexel.edu | Expiration: 12/31/25
================================================================

Setup Species
  * label = C
  * stat_segment_length (b) = 0.14
  * smearing_length (a) = 1.00
  * smearing_length_units = Rg
  * charge = 0.000000
Setup Species
  * label = W
  * stat_segment_length (b) = 1.00
  * smearing_length (a) = 1.00
  * smearing_length_units = Rg
  * charge = 0.000000
Initialize Cell
  * dim = 3
  * length unit = Rg
  * cell_lengths = [99.009901 20.099751 20.099751 ]
  * tilt_factors = [0.000000 0.000000 0.000000 ]
  * cell_tensor:
      [ 99.01,  0.00,  0.00 ]
      [  0.00, 20.10,  0.00 ]
      [  0.00,  0.00, 20.10 ]
Initialize FieldLayout
  * npw = [100 20 20 ]
  * random_seed = 777
  * backend = NOT_SPECIFIED
  * backend = CUDA (chosen by fieldlib)
  * Execution device: [GPU 0] NVIDIA GeForce RTX 2080 Ti (ComputeCapability 7.5)
Setup MoleculePolymerLinear
  * chain_type = discrete
  * label = mol0
Setup MoleculeBead
  * label = mol1
  * volume_frac = 0.50
  * bead_species = W
Setup ModelChiMultiSpecies
  * bref = 1.00
  * Nref = 2.00
  * C = 5.00
  * rho0 = 146.97
  * compressibility_type = Helfand
  * inverse_zetaN = 0.01
  * zeta = 100.00
  * chiN_matrix =
[ [ 0.00, 4.92 ],
  [ 4.92, 0.00 ]]
  * X_matrix =
[ [ 14.14, 14.49 ],
  [ 14.49, 14.14 ]]
  * O_matrix (cols are eigenvectors) =
[ [ -0.71,  0.71 ],
  [  0.71,  0.71 ]]
  * eigenvalues = -0.35 28.63
  * A_matrix =
[ [ -0.71,  0.71 ],
  [  0.71,  0.71 ]]
  * nfields = 2
  * has_electrostatic_field = 0
Initialize MoleculePolymerLinear
  * Linear polymer has single chain:
    - type = discrete
    - nbeads = 16
    - discrete_type = Gaussian
    - sequence_type = monomers
    - nblocks = 1
    - block_fractions =  1.000000
    - block_species =  C
    - monomer_sequence (compact) =  C_16
    - monomer_sequence (long) = CCCCCCCCCCCCCCCC
    - species_fractions =  1.00 0.00
  * volume_frac = 0.500000000000000
  * ncopies = 12500.000000000000000
Initialize MoleculeBead
  * volume_frac = 0.500000000000000
  * ncopies = 200000.000000000000000
Initialize ModelChiMultiSpecies
  * total net_charge / (CV) = 0.000000000000000
  * exchange_field_names = 'mu_1' 'mu_2'
  * initfields = yes
    - reading 'exchange' fields from file 'fields.in'
Initialize DriverSCFTCL
  * is_complex_langevin = false
  * dt = 1.000000000000000
  * output_freq = 1000
  * output_fields_freq = 1000
  * block_size = 1000
  * field_updater = EMPEC
    - predictor_corrector = true
    - update_order = simultaneous
    - lambdas = 0.020000 0.100000
    - adaptive_timestep = false

  * field_stop_tolerance = 0.000010000000000
Calculating Operators:
  * OperatorHamiltonian
    - type = scalar
    - averaging_type = none
    - save_history = true
    - units = beta_Rg3_over_V
Output to files:
  * scalar_operators  -> "operators.dat"
  * species_densities -> "density.dat"
                      ->  save_history = false
  * exchange_fields   -> "fields.dat"
                      ->  save_history = false
  * field_errors      -> "error.dat"

Running DriverSCFTCL
  * nsteps = 100000
  * resume = false
  * store_operator_timeseries = true

Starting Run.
Step:      1  H: -2.505744e+00 -3.209e-16i  FieldErr: 5.43e+00  TPS: 112.69
Step:   1000  H: -2.389431e+00 +2.488e-16i  FieldErr: 3.91e-04  TPS: 143.06
Step:   2000  H: -2.389422e+00 +9.572e-17i  FieldErr: 6.29e-05  TPS: 143.87
Step:   3000  H: -2.389422e+00 -1.444e-16i  FieldErr: 1.41e-05  TPS: 144.61
SCFT converged in 3240 steps to tolerance: 0.00 (fields)
* H = -2.38942 -7.28543e-16i
Run completed in 23.28 seconds (139.12 steps / sec).
-2.389421769785

The units of H1 and H2 are $β H R_{g, r e f}^{3} / V$ where $β = 1 / (k_{B} T)$ and $V = L_{x} A_{y z}$ . We first convert $\frac{d F}{d A}$ to J/m $^{2}$ .

[9]:

beta_dF_over_V = H2 - H1
print(f'{beta_dF_over_V = :.10f}')
kB = 1.38e-23 # J/K
T = 298.15 #K

dF = beta_dF_over_V * Lx * Ayz * kB * T # J
nm2_to_m2 = (1e-9)**2
dF_dA = dF / (dA  * nm2_to_m2) # J / m2
print(f"{dF_dA = } J / m2")

beta_dF_over_V = 0.0020423529
dF_dA = 0.08403199769377069 J / m2

Finally, we can compute $γ$ in units of mJ/m $^{2}$ :

[10]:

J_to_mJ = 1e+3
gamma = 0.5 * dF_dA * J_to_mJ
print(f"{gamma = } mJ / m2")

gamma = 42.015998846885346 mJ / m2

This value of the surface tension is in reasonable agreement with 55.2 mJ/m $^{2}$ as reported in https://pubs.acs.org/doi/full/10.1021/la960800g.