{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b364a736",
   "metadata": {},
   "source": [
    "# Calculating surface tension \n",
    "\n",
    "In this tutorial we demonstrate how to calculate the surface tension between water and hexadecane.\n",
    "\n",
    "The first place to start is the definition of the surface tension $\\gamma$\n",
    "\n",
    "$$ \\frac{dF}{dA} \\equiv 2 \\gamma  $$\n",
    "\n",
    "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. \n",
    "\n",
    "The Helmholtz free energy can be conveniently found using SCFT by locating the saddle point fields $\\mathbf{w}^*$:\n",
    "\n",
    "$$ \n",
    "\\beta F = - \\ln \\mathcal{Z}_c \\approx -\\ln (e^{-H[\\mathbf{w}^*]}) = H[\\mathbf{w}^*].\n",
    "$$\n",
    "\n",
    "To begin, lets consider the units in OpenFTS. OpenFTS provides two options to specify units. \n",
    "\n",
    "1. `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)\n",
    "2. `b` : units are given in terms of a reference statistical segement length $b_{ref}$\n",
    "\n",
    "We will use `Rg` units and will use the convention that $R_{g,ref} = 1$ nanometer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3d53e96b",
   "metadata": {},
   "outputs": [],
   "source": [
    "units = 'Rg'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8636362",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "69d947ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 16 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de6e791",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "18655eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_C = 0.137 # units of bref = 2.45 nm\n",
    "b_W = 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "018096bf",
   "metadata": {},
   "source": [
    "For simplicity, we will also set the smearing length for hexadecane and water to $a_C = a_W = 1$ nm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bf4c3776",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_C = a_W = 1 # nm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4f065ff",
   "metadata": {},
   "source": [
    "Finally, we need to specify the interaction parameters between the C and W monomers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d458a0c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi = 2.460621 # this was found using COSMO-RS\n",
    "zeta = 100 # Arbitrarily specified"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bc4175b",
   "metadata": {},
   "source": [
    "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):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eea54ed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openfts\n",
    "import numpy as np\n",
    "\n",
    "def run(Lx, Ayz, dA=0):\n",
    "\n",
    "    # perturb yz area while keeping total volume constant\n",
    "    Ly = Lz = np.sqrt(Ayz + dA)\n",
    "    Lx = Lx*Ayz/(Ly*Lz)\n",
    "\n",
    "    fts = openfts.OpenFTS()\n",
    "\n",
    "    fts.cell = openfts.Cell(cell_scale=1.0, cell_lengths=[Lx,Ly,Lz], dimension=3,\n",
    "                            tilt_factors=[0, 0, 0], length_unit=units)\n",
    "    fts.field_layout = openfts.FieldLayout(npw=[100, 20, 20], random_seed=777)\n",
    "\n",
    "    scft = openfts.driver.SCFT(dt=1.0, nsteps=100000, output_freq=1000, stop_tolerance=1e-5)\n",
    "    scft.field_updater = openfts.field_updater.EMPEC(update_order='simultaneous',\n",
    "                                                     adaptive_timestep=False, lambdas=[0.02, 0.1])\n",
    "    fts.driver = scft\n",
    "\n",
    "    chiMultiSpecies = openfts.model.ChiMultiSpecies(Nref=2, bref=1.0, chi_array=[chi], zeta=zeta, C=5)\n",
    "    chiMultiSpecies.init_fields = openfts.init_field.FromFile('fields.in', field_type='exchange')\n",
    "    fts.model = chiMultiSpecies\n",
    "\n",
    "    fts.molecules.append(openfts.molecule.PolymerLinearDiscrete(monomer_sequence='C'*N, nbeads=N, volume_frac=0.5))\n",
    "    fts.molecules.append(openfts.molecule.Bead(bead_species='W', volume_frac=0.5))\n",
    "    fts.operators.append(openfts.operator.Hamiltonian(units='beta_Rg3_over_V'))\n",
    "    fts.output = openfts.Output()\n",
    "    fts.species.append(openfts.Species(label='C', stat_segment_length=b_C,\n",
    "                       smearing_length=a_C, smearing_length_units=units))\n",
    "    fts.species.append(openfts.Species(label='W', stat_segment_length=b_W,\n",
    "                       smearing_length=a_W, smearing_length_units=units))\n",
    "    fts.run()\n",
    "\n",
    "    H = openfts.get_scalar_operator('operators.dat', 'H.real', 'final')\n",
    "    return H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf3635f2",
   "metadata": {},
   "source": [
    "Now lets run SCFT for a specified `Lx`, `Ayz` and `dA=0` and save the result in `H1`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1cac12b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================\n",
      "                           OpenFTS                              \n",
      "\n",
      "OpenFTS version: heads/develop f48f964 DIRTY\n",
      "FieldLib version: heads/master c3fc5fd CLEAN\n",
      "FieldLib precision: double\n",
      "Licensee: Joshua Lequieu | jl4347@drexel.edu | Expiration: 12/31/25\n",
      "================================================================\n",
      "\n",
      "Setup Species\n",
      "  * label = C\n",
      "  * stat_segment_length (b) = 0.137\n",
      "  * smearing_length (a) = 1\n",
      "  * smearing_length_units = Rg\n",
      "  * charge = 0.000000\n",
      "Setup Species\n",
      "  * label = W\n",
      "  * stat_segment_length (b) = 1\n",
      "  * smearing_length (a) = 1\n",
      "  * smearing_length_units = Rg\n",
      "  * charge = 0.000000\n",
      "Initialize Cell\n",
      "  * dim = 3\n",
      "  * length unit = Rg\n",
      "  * cell_lengths = [100.000000 20.000000 20.000000 ]\n",
      "  * tilt_factors = [0.000000 0.000000 0.000000 ]\n",
      "  * cell_tensor: \n",
      "      [ 100,   0,   0 ]\n",
      "      [   0,  20,   0 ]\n",
      "      [   0,   0,  20 ]\n",
      "Initialize FieldLayout\n",
      "  * npw = [100 20 20 ]\n",
      "  * random_seed = 777\n",
      "  * backend = NOT_SPECIFIED\n",
      "  * backend = CUDA (chosen by fieldlib)\n",
      "  * Execution device: [GPU 0] NVIDIA GeForce RTX 2080 Ti (ComputeCapability 7.5)\n",
      "Setup MoleculePolymerLinear\n",
      "  * chain_type = discrete\n",
      "  * label = mol0\n",
      "Setup MoleculeBead\n",
      "  * label = mol1\n",
      "  * volume_frac = 0\n",
      "  * bead_species = W\n",
      "Setup ModelChiMultiSpecies\n",
      "  * bref = 1\n",
      "  * Nref = 2\n",
      "  * C = 5\n",
      "  * rho0 = 146.969\n",
      "  * compressibility_type = Helfand\n",
      "  * inverse_zetaN = 0.005\n",
      "  * zeta = 100\n",
      "  * chiN_matrix = \n",
      "[ [       0, 4.92124 ],\n",
      "  [ 4.92124,       0 ]]\n",
      "  * X_matrix = \n",
      "[ [ 14.1421, 14.4901 ],\n",
      "  [ 14.4901, 14.1421 ]]\n",
      "  * O_matrix (cols are eigenvectors) = \n",
      "[ [ -0.707107,  0.707107 ],\n",
      "  [  0.707107,  0.707107 ]]\n",
      "  * eigenvalues = -0.347984 28.6323 \n",
      "  * A_matrix = \n",
      "[ [ -0.707107,  0.707107 ],\n",
      "  [  0.707107,  0.707107 ]]\n",
      "  * nfields = 2\n",
      "  * has_electrostatic_field = 0\n",
      "Initialize MoleculePolymerLinear\n",
      "  * Linear polymer has single chain: \n",
      "    - type = discrete\n",
      "    - nbeads = 16\n",
      "    - discrete_type = Gaussian\n",
      "    - sequence_type = monomers\n",
      "    - nblocks = 1\n",
      "    - block_fractions =  1.000000\n",
      "    - block_species =  C\n",
      "    - monomer_sequence (compact) =  C_16\n",
      "    - monomer_sequence (long) = CCCCCCCCCCCCCCCC\n",
      "    - species_fractions =  1 0\n",
      "  * volume_frac = 0.5\n",
      "  * ncopies = 12500\n",
      "Initialize MoleculeBead\n",
      "  * volume_frac = 0.5\n",
      "  * ncopies = 200000\n",
      "Initialize ModelChiMultiSpecies\n",
      "  * total net_charge / (CV) = 0\n",
      "  * exchange_field_names = 'mu_1' 'mu_2' \n",
      "  * initfields = yes\n",
      "    - reading 'exchange' fields from file 'fields.in'\n",
      "Initialize DriverSCFTCL\n",
      "  * is_complex_langevin = false\n",
      "  * dt = 1\n",
      "  * output_freq = 1000\n",
      "  * output_fields_freq = 1000\n",
      "  * block_size = 1000\n",
      "  * field_updater = EMPEC\n",
      "    - predictor_corrector = true\n",
      "    - update_order = simultaneous\n",
      "    - lambdas = 0.020000 0.100000 \n",
      "    - adaptive_timestep = false\n",
      "\n",
      "  * field_stop_tolerance = 1e-05\n",
      "Calculating Operators:\n",
      "  * OperatorHamiltonian\n",
      "    - type = scalar\n",
      "    - averaging_type = none\n",
      "    - save_history = true\n",
      "    - units = beta_Rg3_over_V\n",
      "Output to files: \n",
      "  * scalar_operators  -> \"operators.dat\"\n",
      "  * species_densities -> \"density.dat\"\n",
      "                      ->  save_history = false\n",
      "  * exchange_fields   -> \"fields.dat\"\n",
      "                      ->  save_history = false\n",
      "  * field_errors      -> \"error.dat\"\n",
      "\n",
      "Running DriverSCFTCL\n",
      "  * nsteps = 100000\n",
      "  * resume = false\n",
      "  * store_operator_timeseries = true\n",
      "\n",
      "Starting Run.\n",
      "Step:      1  H: -2.509370e+00 -2.609e-16i  FieldErr: 5.44e+00  TPS:  59.58\n",
      "Step:   1000  H: -2.391474e+00 +1.748e-17i  FieldErr: 3.87e-04  TPS: 141.68\n",
      "Step:   2000  H: -2.391465e+00 -1.947e-16i  FieldErr: 6.19e-05  TPS: 144.59\n",
      "Step:   3000  H: -2.391464e+00 +5.016e-16i  FieldErr: 1.39e-05  TPS: 142.83\n",
      "SCFT converged in 3230 steps to tolerance: 0.00 (fields)\n",
      "* H = -2.39146 -2.20678e-16i\n",
      "Run completed in 23.20 seconds (139.18 steps / sec).\n",
      "-2.391464122693\n"
     ]
    }
   ],
   "source": [
    "Lx = 100 #nm \n",
    "Ayz = 400 #nm\n",
    "H1 = run(Lx, Ayz, dA=0.0)\n",
    "print(H1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c964f0c",
   "metadata": {},
   "source": [
    "Next lets run SCFT for a specified `Lx`, `Ayz` and `dA=4` and save the result in `H2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ea290067",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================\n",
      "                           OpenFTS                              \n",
      "\n",
      "OpenFTS version: heads/develop f48f964 DIRTY\n",
      "FieldLib version: heads/master c3fc5fd CLEAN\n",
      "FieldLib precision: double\n",
      "Licensee: Joshua Lequieu | jl4347@drexel.edu | Expiration: 12/31/25\n",
      "================================================================\n",
      "\n",
      "Setup Species\n",
      "  * label = C\n",
      "  * stat_segment_length (b) = 0.14\n",
      "  * smearing_length (a) = 1.00\n",
      "  * smearing_length_units = Rg\n",
      "  * charge = 0.000000\n",
      "Setup Species\n",
      "  * label = W\n",
      "  * stat_segment_length (b) = 1.00\n",
      "  * smearing_length (a) = 1.00\n",
      "  * smearing_length_units = Rg\n",
      "  * charge = 0.000000\n",
      "Initialize Cell\n",
      "  * dim = 3\n",
      "  * length unit = Rg\n",
      "  * cell_lengths = [99.009901 20.099751 20.099751 ]\n",
      "  * tilt_factors = [0.000000 0.000000 0.000000 ]\n",
      "  * cell_tensor: \n",
      "      [ 99.01,  0.00,  0.00 ]\n",
      "      [  0.00, 20.10,  0.00 ]\n",
      "      [  0.00,  0.00, 20.10 ]\n",
      "Initialize FieldLayout\n",
      "  * npw = [100 20 20 ]\n",
      "  * random_seed = 777\n",
      "  * backend = NOT_SPECIFIED\n",
      "  * backend = CUDA (chosen by fieldlib)\n",
      "  * Execution device: [GPU 0] NVIDIA GeForce RTX 2080 Ti (ComputeCapability 7.5)\n",
      "Setup MoleculePolymerLinear\n",
      "  * chain_type = discrete\n",
      "  * label = mol0\n",
      "Setup MoleculeBead\n",
      "  * label = mol1\n",
      "  * volume_frac = 0.50\n",
      "  * bead_species = W\n",
      "Setup ModelChiMultiSpecies\n",
      "  * bref = 1.00\n",
      "  * Nref = 2.00\n",
      "  * C = 5.00\n",
      "  * rho0 = 146.97\n",
      "  * compressibility_type = Helfand\n",
      "  * inverse_zetaN = 0.01\n",
      "  * zeta = 100.00\n",
      "  * chiN_matrix = \n",
      "[ [ 0.00, 4.92 ],\n",
      "  [ 4.92, 0.00 ]]\n",
      "  * X_matrix = \n",
      "[ [ 14.14, 14.49 ],\n",
      "  [ 14.49, 14.14 ]]\n",
      "  * O_matrix (cols are eigenvectors) = \n",
      "[ [ -0.71,  0.71 ],\n",
      "  [  0.71,  0.71 ]]\n",
      "  * eigenvalues = -0.35 28.63 \n",
      "  * A_matrix = \n",
      "[ [ -0.71,  0.71 ],\n",
      "  [  0.71,  0.71 ]]\n",
      "  * nfields = 2\n",
      "  * has_electrostatic_field = 0\n",
      "Initialize MoleculePolymerLinear\n",
      "  * Linear polymer has single chain: \n",
      "    - type = discrete\n",
      "    - nbeads = 16\n",
      "    - discrete_type = Gaussian\n",
      "    - sequence_type = monomers\n",
      "    - nblocks = 1\n",
      "    - block_fractions =  1.000000\n",
      "    - block_species =  C\n",
      "    - monomer_sequence (compact) =  C_16\n",
      "    - monomer_sequence (long) = CCCCCCCCCCCCCCCC\n",
      "    - species_fractions =  1.00 0.00\n",
      "  * volume_frac = 0.500000000000000\n",
      "  * ncopies = 12500.000000000000000\n",
      "Initialize MoleculeBead\n",
      "  * volume_frac = 0.500000000000000\n",
      "  * ncopies = 200000.000000000000000\n",
      "Initialize ModelChiMultiSpecies\n",
      "  * total net_charge / (CV) = 0.000000000000000\n",
      "  * exchange_field_names = 'mu_1' 'mu_2' \n",
      "  * initfields = yes\n",
      "    - reading 'exchange' fields from file 'fields.in'\n",
      "Initialize DriverSCFTCL\n",
      "  * is_complex_langevin = false\n",
      "  * dt = 1.000000000000000\n",
      "  * output_freq = 1000\n",
      "  * output_fields_freq = 1000\n",
      "  * block_size = 1000\n",
      "  * field_updater = EMPEC\n",
      "    - predictor_corrector = true\n",
      "    - update_order = simultaneous\n",
      "    - lambdas = 0.020000 0.100000 \n",
      "    - adaptive_timestep = false\n",
      "\n",
      "  * field_stop_tolerance = 0.000010000000000\n",
      "Calculating Operators:\n",
      "  * OperatorHamiltonian\n",
      "    - type = scalar\n",
      "    - averaging_type = none\n",
      "    - save_history = true\n",
      "    - units = beta_Rg3_over_V\n",
      "Output to files: \n",
      "  * scalar_operators  -> \"operators.dat\"\n",
      "  * species_densities -> \"density.dat\"\n",
      "                      ->  save_history = false\n",
      "  * exchange_fields   -> \"fields.dat\"\n",
      "                      ->  save_history = false\n",
      "  * field_errors      -> \"error.dat\"\n",
      "\n",
      "Running DriverSCFTCL\n",
      "  * nsteps = 100000\n",
      "  * resume = false\n",
      "  * store_operator_timeseries = true\n",
      "\n",
      "Starting Run.\n",
      "Step:      1  H: -2.505744e+00 -3.209e-16i  FieldErr: 5.43e+00  TPS: 112.69\n",
      "Step:   1000  H: -2.389431e+00 +2.488e-16i  FieldErr: 3.91e-04  TPS: 143.06\n",
      "Step:   2000  H: -2.389422e+00 +9.572e-17i  FieldErr: 6.29e-05  TPS: 143.87\n",
      "Step:   3000  H: -2.389422e+00 -1.444e-16i  FieldErr: 1.41e-05  TPS: 144.61\n",
      "SCFT converged in 3240 steps to tolerance: 0.00 (fields)\n",
      "* H = -2.38942 -7.28543e-16i\n",
      "Run completed in 23.28 seconds (139.12 steps / sec).\n",
      "-2.389421769785\n"
     ]
    }
   ],
   "source": [
    "dA = 4 # nm2\n",
    "H2 = run(Lx, Ayz, dA=dA)\n",
    "print(H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d798dd9",
   "metadata": {},
   "source": [
    "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$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "880dda3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_dF_over_V = 0.0020423529\n",
      "dF_dA = 0.08403199769377069 J / m2\n"
     ]
    }
   ],
   "source": [
    "beta_dF_over_V = H2 - H1\n",
    "print(f'{beta_dF_over_V = :.10f}')\n",
    "kB = 1.38e-23 # J/K\n",
    "T = 298.15 #K\n",
    "\n",
    "dF = beta_dF_over_V * Lx * Ayz * kB * T # J\n",
    "nm2_to_m2 = (1e-9)**2\n",
    "dF_dA = dF / (dA  * nm2_to_m2) # J / m2\n",
    "print(f\"{dF_dA = } J / m2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab72833d",
   "metadata": {},
   "source": [
    "Finally, we can compute $\\gamma$ in units of mJ/m$^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9eb2278e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma = 42.015998846885346 mJ / m2\n"
     ]
    }
   ],
   "source": [
    "J_to_mJ = 1e+3\n",
    "gamma = 0.5 * dF_dA * J_to_mJ\n",
    "print(f\"{gamma = } mJ / m2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f05e4e1",
   "metadata": {},
   "source": [
    "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`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (base env)",
   "language": "python",
   "name": "base"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}