Grand Canonical Monte Carlo of atomic species using tabulated pair potentials

In [2]:

    

from __future__ import division, unicode_literals, print_function
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np, pandas as pd
import os.path, os, sys, json, filecmp, copy
plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})
try:
    workdir
except NameError:
    workdir=%pwd
else:
    %cd $workdir
print(workdir)


    

    




/Users/mikael/github/faunus-notebooks/tabulated-potential

In [7]:

    

%%bash -s "$workdir"
cd $1

if [ ! -d "faunus/" ]; then
    git clone https://github.com/mlund/faunus.git
    cd faunus
    #git checkout a0f0b46
else
    cd faunus
fi

# if different, copy custom grand.cpp into faunus
if ! cmp ../grand.cpp src/examples/grand.cpp >/dev/null 2>&1
then
    cp ../grand.cpp src/examples/
fi

pwd
CXX=clang++ CC=clang cmake . -DCMAKE_BUILD_TYPE=RelWithDebInfo -DENABLE_APPROXMATH=on &>/dev/null
make example_grand -j4


    

    




/Users/mikael/github/faunus-notebooks/tabulated-potential/faunus
[ 36%] Built target xdrfile
[ 81%] Built target libfaunus
Scanning dependencies of target example_grand
[ 90%] Building CXX object src/examples/CMakeFiles/example_grand.dir/grand.cpp.o
[100%] Linking CXX executable grand
[100%] Built target example_grand

In [4]:

    

r = np.arange(0, 200, 0.1)
u = np.sin(r)*np.exp(-r/6)
plt.plot(r, u)
d = np.dstack((r,u))
np.savetxt('cation-anion.dat', np.c_[r,u])


    

In [10]:

    

%cd $workdir

def mkinput():
    js = {
        "atomlist" : {
            "cat" : { "eps": 0.15, "sigma":4.0, "dp":40, "activity":Cs, "q":1.0 },
            "an"  : { "eps": 0.20, "sigma":4.0, "dp":10, "activity":Cs, "q":-1.0 }
        },
        "moleculelist" : { 
            "salt" : { "atoms":"cat an", "atomic":True, "Ninit":50 }
        },
        "moves" : {
            "atomtranslate" : {
                "salt" : { "peratom":True }
            },
            "atomgc": { "molecule": "salt" }
        },
        "energy" : {
            "nonbonded" : {
                "pairpotentialmap" : {
                    "spline" : { "rmin":4.0, "rmax":100, "utol":0.01  },
                    "cat an" : {
                        "fromdisk" : "../cation-anion.dat",
                        "_coulomb" : { "epsr": 80.0 }
                    },
                    "default" : {
                        "lennardjones" : {}
                    }
                }
            }
        },
        "system" : {
            "temperature" : 298.15,
            "cuboid" : { "len":50 },
            "mcloop": { "macro": 10, "micro": micro }
        }
    }
    with open('gc.json', 'w+') as f:
        f.write(json.dumps(js, indent=4))

Cs_range = [0.1]
        
for Cs in Cs_range:
    pfx='Cs'+str(Cs)
    if True: #not os.path.isdir(pfx):
        %mkdir -p $pfx
        %cd $pfx
        # equilibration run (no translation)
        !rm -fR state
        micro=1000
        mkinput()
        !../faunus/src/examples/grand > eq
        
        # production run
        micro=1000
        mkinput()
        %time !../faunus/src/examples/grand
        %cd ..
print('done.')


    

    




/Users/mikael/github/faunus-notebooks/tabulated-potential
/Users/mikael/github/faunus-notebooks/tabulated-potential/Cs0.1
State file not found.
Warning: energy change from move returns not-a-number (NaN)
Warning: energy change from move returns not-a-number (NaN)
Reading space state file 'state'. OK!
  Resizing particle vector from 100 --> 16.
  Read 16 particle(s).
  Read 1 group(s).
  Restoring random number generator state.

 .................
  Atom Properties  
 *****************
  Number of entries:       3
  Element info:
    unk   an    cat   

 ...............................
  Simulation Space and Geometry  
 *******************************
  Boundary                 Cuboid
  Volume                   125000 ų = 125 nm³ = 1.25e-22 liters
  Sidelengths              50 50 50 (Å)
  Scale directions         XYZ
  Number of particles      16
  Electroneutrality        Yes 0
  System sanity check      Passed
  Number of molecule types 1
  Groups:
    1     [0-15]           salt        N/V = 0.000128 Å⁻³ = 212.549 mM  

 ..........................................
  Energy: Nonbonded N² - pairpotentialmap  
 ******************************************
  Rmin                4
  Rmax                100
  Utol                0.01
  Ftol                -1
  Umaxtol             -1
  Fmaxtol             -1

Lennard-Jones w. Lorentz-Berthelot Mixing:
an-an cat-cat 

fromdisk:  [0:199.9]:
an-cat 


  Steps:    1 / 1000        Macrosteps/min: 3750    ETA: Tue Mar 15 16:14:23 2016
  Steps:    2 / 2000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    3 / 3000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    4 / 4000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    5 / 5000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    6 / 6000        Macrosteps/min: 3333    ETA: Tue Mar 15 16:14:23 2016
  Steps:    7 / 7000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    8 / 8000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:    9 / 9000        Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
  Steps:   10 / 10000       Macrosteps/min: 3529    ETA: Tue Mar 15 16:14:23 2016
Writing to file 'confout.pqr'. OK!
  Writing space state file 'state'. OK!

 ...................
  MC Steps and Time  
 *******************
  Steps (macro micro tot)  10∙1000 = 10000
  Time elapsed             0 min = 0 h

 .........................
  System Energy and Drift  
 *************************
  Average                  -0.290931 kT, σ=0.39317
  Initial energy           -0.501547 kT
  Initial + changes        -0.316172 kT
  Total energy drift       -4.219e-15 kT (1.334e-12﹪)

 ....................................
  Markov Move: P R O P A G A T O R S  
 ************************************

 ...................................
  Markov Move: Grand Canonical Salt  
 ***********************************
  Number of trials              5021
  Relative time consumption     0.151814
  Acceptance                    76.3394﹪
  Runfraction                   100﹪
  Total energy change           -22.3806 kT
  Number of GC species          

    Ion       activity  ⟨c/M⟩     ⟨ɣ±⟩      ⟨N⟩       
    an        0.1       0.092927  1.0761    6.9952    
    cat       0.1       0.092927  1.0761    6.9952    

 ..........................................
  Markov Move: Single Particle Translation  
 ******************************************
  Number of trials              69500
  Relative time consumption     0.715403
  Acceptance                    90.6432﹪
  Runfraction                   100﹪
  Total energy change           22.566 kT
  Average moves/particle        4633
  Displacement vector           1 1 1

  Individual particle movement:

           dp    Acc. ﹪     Nmoves      ⟨r²⟩/Ų     √⟨r²⟩/Å
    an     10    91.3        34654       22.7        4.77        
    cat    40    90          34846       360         19          

 .........................................
  Analysis: Multi Particle Widom Analysis  
 *****************************************
  Reference:                    doi:10/dkv4s6
  Sample interval               0
  Number of sample events       10000
  Number of insertions          10000
  Excess chemical pot.          -0.000718813 kT
  Mean activity coefficient     0.999281
  Ghost particles               an cat 
CPU times: user 4.28 ms, sys: 5.56 ms, total: 9.84 ms
Wall time: 312 ms
/Users/mikael/github/faunus-notebooks/tabulated-potential
done.

In [ ]: