I am trying to define a region as a slab and add force to solvent molecules found in this slab in order to create a flow. In this way, I intend to examine the solvent transport through a channel via non-equilibrium MD. I tried this:
group gSOLVENT type 5:14
region top block EDGE EDGE EDGE EDGE EDGE -25.8 units box side in
fix AppForce1 gSOLVENT addforce 0.0 0.0 0.1 region top
It obviously works but some solvent molecules straddle the region boundary and I want to apply force to all atoms of a molecule (even if only one atom of the molecule falls within the slab). For this, I need to define a dynamic group including both my original in-slab atoms and other atoms having the same molecule ID with in-slab atoms. I am aware of that “include molecule” keyword in group command does not work in a dynamic way. After searching the manual and previous topics, I tried using a combination of a “compute” and an “atom-style variable” as follows:
group gSOLVENT type 5:14
region top block EDGE EDGE EDGE EDGE EDGE -25.8 units box side in
group gSOLVENTTOP dynamic gSOLVENT region top
compute cTOPID gSOLVENTTOP property/atom mol # produces a per-atom vector containing molecule ids of my in-slab atoms
variable vTOPMOL atom “mol==c_cTOPID” # expected to define another per-atom vector of zeros and ones, defining “one” for all atoms having molecule IDs of interest
It didn’t add the atoms I wanted to add. What do I miss and what else can I try? I am using the version 29 Sep 2021. I am also adding my input and data files:
I am not sure how physical your scheme is – but if you’re just simulating water, why not just apply force to the water oxygens and let SHAKE worry about distributing that force to the hydrogens?
Dear @srtee, thank you so much for your answer. I am also simulating organic solvents, larger molecules like hexane and heptane as well as water. First, I started with water using the method you suggested and there was no problem. However, when I switched to other solvents, the force that I applied became insufficient. We suspected that problem may be related to nonuniform force distribution within the solvent molecules. As molecule gets larger, only a part of it falls within the slab.
Please also factor in Newton’s second law of motion. If you want to see the same acceleration for different size objects with different (total) masses, you need proportional forces. When you apply a constant size force to multiple atoms the total force is the sum of the individual forces.
Yes, which is why this “slab” approach does not seem physically trustworthy to me. There are other established approaches in literature to do non-equilibrium molecular dynamics.
You could try to set up a per-atom variable which is 1 in the slab region, 0 far outside, and gradually changes from 1 to 0 in between in some transition region, and then specify that the added force (or acceleration, using fix gravity) is your desired magnitude times that variable. That might help if you’re concerned about abruptly changing forces.
Dear @akohlmey, thank you so much for your warning. Actually, I don’t want to see the same acceleration for different solvent molecules (e.g. water and hexane). I want to apply the same pressure difference and observe different transport behaviours. In NEMD studies that I follow, pressure difference is calculated as deltaP=n.f/A, n is the number of particles, f is the force and A is the channel area. So I will use an f constant that is proportional to “n”, trying to even out n.f for different solvent systems. As a result, I expect different acceleration under the same pressure difference. My worry is to introduce an additional, nonphysical difference originating from the molecules that are partially outside the slab. This problem may not be crucial for molecules like water but we thought its effect may be pronounced for large molecules.
Dear @srtee, did you mean using a piston by “established approaches”? This is what we recently thought as an alternative. Also, many thanks for your suggestion, we will try the transition region.
But that is exactly what I am saying. If you apply the force to the “central” atom of a molecule instead of all its atoms, you have to multiply it with the number of atoms per molecule to have the same total effect.
For the sake of completeness and to give an instructive example, here is some LAMMPS input fragment, that can be used to select all atoms of a molecule, if only one of them is contained in a region (or conforms to some other condition coded in the insel atom style variable).
region sel sphere 0.0 0.0 0.0 10.0 units box # region of interest
variable insel atom rmask(sel) # flag atoms in region of interest with 1, others 0
compute cpa all chunk/atom molecule # define per-molecule chunks
compute mpc all reduce/chunk cpa max v_insel # all chunks with 1+ atom in region are 1, others 0
compute spm all chunk/spread/atom cpa c_mpc # spread per-chunk value to all atoms in chunk
After these steps a reference to c_spm can be used for further processing, e.g. to be referenced in another atom style variable to define a dynamic group.
Dear @akohlmey, we are so grateful for your help. Thanks to your input fragment, now I am able to do what I exactly asked for. I made a small addition and defined the force as a variable to make “n.f” constant. I am sharing below the last version of the related part in case anyone needs it in future.
group gSOLVENT type 5:14
region top block EDGE EDGE EDGE EDGE EDGE -25.8 units box side in
variable insel atom rmask(top) # flag atoms in region of interest with 1, others 0
compute cpa all chunk/atom molecule # define per-molecule chunks
compute mpc all reduce/chunk cpa max v_insel # all chunks with 1+ atom in region are 1, others 0
compute spm all chunk/spread/atom cpa c_mpc # spread per-chunk value to all atoms in chunk
variable vspm atom c_spm # reference the compute in another atom style variable
group gSLABTOPMOL dynamic gSOLVENT var vspm # define the dynamic group
variable numslabtmol equal count(gSLABTOPMOL) # number of atoms in-slab + include molecule