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 \(\gamma\)
\[\frac{dF}{dA} \equiv 2 \gamma\]
where \(F\) is the Helmholtz free energy and \(A\) is the area of the interface. Our approach will be to calculate the \(dF/dA\) derivative using a finite difference.
The Helmholtz free energy can be conveniently found using SCFT by locating the saddle point fields \(\mathbf{w}^*\):
\[\beta F = - \ln \mathcal{Z}_c \approx -\ln (e^{-H[\mathbf{w}^*]}) = H[\mathbf{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,ref}= b_{ref}\sqrt{N_{ref}/6}\) (continuous chains) or \(R_{g,ref} = b_{ref}\sqrt{(N_{ref}-1)/6}\) (discrete chains)b: units are given in terms of a reference statistical segement length \(b_{ref}\)
We will use Rg units and will use the convention that \(R_{g,ref} = 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_{yz}\) (yz cross sectional area) and \(dA\) (perturbation in cross sectional area):
[6]:
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(cell_scale=1.0, cell_lengths=[Lx,Ly,Lz], dimension=3,
tilt_factors=[0, 0, 0], length_unit=units)
fts.field_layout(npw=[100, 20, 20], random_seed=777)
fts.driver(dt=1.0, nsteps=100000, output_freq=1000, stop_tolerance=1e-5,
type='SCFT')
fts.field_updater(type='EMPEC', update_order='simultaneous',
adaptive_timestep=False, lambdas=[0.02, 0.1])
fts.model(Nref=2, bref=1.0, chi_array=[chi], zeta=zeta, C=5,
type='ChiMultiSpecies')
fts.init_fields_from_file('fields.in', field_type='exchange')
fts.add_molecule(type='PolymerLinear', monomer_sequence='C'*N, nbeads=N,
chain_type='discrete', volume_frac=0.5)
fts.add_molecule(type='Bead', bead_species='W', volume_frac=0.5)
fts.add_operator(type='Hamiltonian', units='beta_Rg3_over_V')
fts.output_default()
fts.add_species(label='C', stat_segment_length=b_C,
smearing_length=a_C, smearing_length_units=units)
fts.add_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 \(\beta H R_{g,ref}^3 / V\) where \(\beta = 1 /(k_B T)\) and \(V = L_x A_{yz}\). We first convert \(\frac{dF}{dA}\) 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 \(\gamma\) 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.