# Free energy of insertion

Thanks Dr. Axel and Carlos what I am trying to do here is the method described in this paper.
Sorry for being vague. I have done the survey. This is the only one I have so far that is closely relevant to the one which I am currently working on. This method is similar to the what carlos had sent involves discretization of the domain…

I am trying to implement this method. They have implement this in lammps. I am not sure how though. They have worked with a very small number of atoms I see now its due to the computational effort it takes to calculate the chemical potential.

The method in a nutshell:

1. Discretize the simulation domain in tiny cubes.
2. Perform insertion in each of the cube, calculate the boltzmannn factor. If its greater than epsilon its classified as a hole.
3. Then in each hole dump 1000 water molecules.
4. Recalcuate average boltzmann factor and calculate excess chemical potential.
5. Repeat this for 10000 frames.
6. Get ensemble average (NPT) of the chemical potential (Expressions given in the paper and Allen and Tildesley)

They are increasing the weight percentage of water slowly until the chemical potential of water in the network is equal to the chemical potential of water outside of the network.

The part where I am stuck is the second step.
All I need is fast insertion in to the boxes after discretization of the domain (which can be done by putting fix deposit in a loop and changing the region of deposition at every iteration). But fix deposit builds neighbour list every step of the deposition process which i think is very very computationally expensive.
That’s why I have asked in the previous post whether fix deposit can acknowledge neigh modify once. But the answer was no.
I got around this problem by using the include command in neigh modify where the “group” is an empty group with which I was able to deposit molecules in to the boxes with very less computational effort.

After getting this configuration
The exclude command in neigh modify can then be used to switch of the interactions in the network, loop over each of the boxes and the find the hole, if its a hole deposit 1000 water molecules (I am not sure how to do this I am guessing using fix pour and fix gravity with a magnitude of zero or fix deposit itself).
Again the exclude command loops over all the atoms in the system and switches off the pair wise interactions in the network
What I rather need is the neighbour list build only for the water molecules. I thought the include command in neigh modify will do that. It doesn’t which was clarified from the previous post.

Any other suggestions to reduce computational effort would be helpful or if what I am trying to do doesn’t make sense let me know

Dr. Axel’s suggestion:

“if you want to find the location available with the lowest energy, you
need to make a huge number of insertion attempts (and then there is
the question if you want to relax the matrix or not. at that point you
could, for example, just as well run a grand canonical ensemble and
adjust the chemical potential until you get a specific desired
occupation.”

Ans: Working on it

Thank you .

jp105210y.pdf (1.17 MB)

I am trying to implement this method. They have implement this in lammps. I am not sure how though.

Why don’t you send an email to the authors? That’s
why the author contact info is included in papers.

Steve

Thanks Dr. Axel and Carlos what I am trying to do here is the method
described in this paper.
Sorry for being vague. I have done the survey. This is the only one I have
so far that is closely relevant to the one which I am currently working on.
This method is similar to the what carlos had sent involves discretization
of the domain.......

I am trying to implement this method. They have implement this in lammps. I

that is not what the paper says. all it says is that they have used
LAMMPS for the MD simulations. it says nothing how the insertions were
performed. most likely with an external software. possibly using the
LAMMPS library interface or just external scripting. the best way to
find out, is to ask the authors.

there is no currently existing feature that does the excluded volume
pre-screening to quickly determine unlikely insertions and thus avoid
the expensive computations for those cases.

axel.

Arun,

You got Steve and Axel’s clear answers. I will go a step further to point to the fact that you can probably use the Voronoi package already in lammps to get a partitioning of the space which “may” be better (with extra help) to the cubic grid-type the authors employed in that paper. Of course this is only and idea and if what you are trying to do is to exactly reproduce the steps of that paper it won’t help. Because lammps can give you the Voronoi volume/atom; before inserting any molecule you could sort the atoms in the system according to their Voronoi volume and a size analysis may help you filter out those where a water molecule would not fit. Then during the insertion step you’ll have to insert the target water mol in the vicinity of each each one of the important atoms following a shell like 3D geometry with some predefined inner and outer radius and compute your exponential factor for the further chemical potential calculation. This way you only need to know the important atom locations not the Voronoi cell structure which the lammps compute does not give you unless your C++ is good enough to go struggle with the Voro++ interface.

Hope you get my point on how to try to tackle this problem only using lammps features.

I cannot fully guarantee the success of this idea so just take it with caution and if have the desire/time try exploring it.

Carlos