Environment similarity with atoms rotating around lattice sites

[Mathieu-Istas]

Dear plumed Users and developers,
I want to study a liquid-solid phase transition, like in this tutorial .
The system I am studying is solid hydrogen, in an HCP phase where each site is occupied by a H2 molecule. I am not sure how to use the environment similarity as a CV, since in simulations H2 molecules are rotating around the HCP lattice sites.
I did try to use a reference pdb file, but this does not result in a CV that can distinguish the liquid and solid phases.
I also tried to use the HCP structure as
cv: ENVIRONMENTSIMILARITY ... SPECIES=1-96 SIGMA=0.023 CRYSTAL_STRUCTURE=HCP LATTICE_CONSTANTS=0.17,0.25 MORE_THAN={CUBIC D_0=0.35 D_MAX=0.56} MEAN

running simulations with OPES_METAD, I did not observe many transitions between liquid and solid. Maybe the overlap is too big, here are the overlaps I obtained for different values of sigma
overlap_plumed-base.dat.pdf
How can I improve the CV taking into account the rotational degree of freedom or using the center of mass of neighboring atoms?

I know we can compute the center of mass in plumed (https://www.plumed.org/doc-v2.9/user-doc/html/_c_o_m.html), however I don't know how to pair atoms ( it is possible to write some code using for instance the Hungarian algorithm, but I don't know how to patch it to plumed.

Thanks,
Mathieu

[PabloPiaggi]

Hi @Mathieu-Istas, I suggest to try using the center of mass of the H2 molecules as the input for the SPECIES keyword. This should work provided that bonds don't break, i.e., you are dealing with a plastic phase and not with a superionic phase. For instance,

c1: COM ATOMS=1-2
c2: COM ATOMS=3-4
.
.
.
c48: COM ATOMS=95-96

cv: ENVIRONMENTSIMILARITY SPECIES=c1,...,c48 SIGMA=0.023 CRYSTAL_STRUCTURE=HCP  LATTICE_CONSTANTS=0.17,0.25 MORE_THAN={CUBIC D_0=0.35 D_MAX=0.56} MEAN

You can find which atoms are forming each molecule using neighbor lists in ASE.

Best,

Pablo

[Mathieu-Istas]

Hi Pablo,
thanks for the answer. I am not sure to understand everything. hydrogen is quite light, and protons often jump from one molecule to another. Is it possible during the LAMMPS simulation to find nearest neighbors (using ASE or something I wrote) and then update the plumed.dat file realtime and to compute the CV? If so is there an example in the documentation that does it? I might have missed that.

Best,
Mathieu

[PabloPiaggi]

I do not think you can define dynamic groups inside plumed. This would be needed to compute the COM of H2 molecules which exchange protons during the simulation. Maybe @gtribello can confirm.

You could consider adding a restraint in LAMMPS to avoid breaking bonds, if this doesn't affect significantly the thermodynamic stability of your solid.

[Mathieu-Istas]

I'm a bit reluctant to fix bonds for the lightest nucleus. I guess I could also use the structure factor as in your silica article. There is a branch there https://github.com/invemichele/plumed2/tree/structure_factor, but it is quite behind the main branch. Do you know if it is possible to patch it to lammps and use it?
Thanks again for the answers

[PabloPiaggi]

Yes, the structure factor in that repo should work well, even if it is quite behind the main branch. I suggest you use it with the LOAD keyword. See for instance this egg: https://www.plumed-nest.org/eggs/21/005/data/CO2/plumed.dat.html

[Mathieu-Istas]

Sorry to bother you again, but I am struggling to use it. If I load the cpp file, I get the error 'lammps/build/plumed_build-prefix/lib/plumed/scripts/mklib.sh ' : No such file or directory. If I try to compile the file with plumed mklib instead, I get the error plumed2/src/structurefactor/StructureFactor_sphericallyAveraged.so: undefined symbol: _ZTIN4PLMD6ActionE. My guess is that this is due to a mismatch between the plumed version used in LAMMPS and the other one. I did build lammps with -D PKG_PLUMED=ON, is there a flag to use a specific plumed version when building lammps?

[gtribello]

Hi all

Could you do something like this for the centers of mass of the hydrogens:

# Calculate distance between first hydrogen and all other 99 hydrogens in the system
d: DISTANCE ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ...
# Apply a switching function on all the distances
s: LESS_THAN ARG=d SWITCH={EXP R_0=0.1}
# Calculate the vectors corresponding to the distances above
dc: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=1,3 ATOMS3=1,4 ATOMS4=1,5 ...
# Multiply the vectors by the switching function - the switching function is 0 for all but the closest atom
vx: CUSTOM ARG=dc.x,s FUNC=x*y PERIODIC=NO
vy: CUSTOM ARG=dc.y,s FUNC=x*y PERIODIC=NO
vz: CUSTOM ARG=dc.z,s FUNC=x*y PERIODIC=NO
# Get the position of the first atom
p1: POSITION ATOM=1
# Now we get the center of mass
comx: CUSTOM ARG=p1.x,vx FUNC=x+0.5*y PERIODIC=NO
comy: CUSTOM ARG=p1.y,vy FUNC=x+0.5*y PERIODIC=NO
comz: CUSTOM ARG=p1.z,vz FUNC=x+0.5*y PERIODIC=NO
# And convert these three arguments to an atom 
m: CONSTANT VALUE=2
q: CONSTANT VALUE=0
com1: ARGS2VATOM XPOS=comx YPOS=comy ZPOS=comz MASS=m CHARGE=q 

If you have got to do this for a lot of atoms you can use CONTACT_MATRIX instead for the distances and the same ideas. You can take a look at the new manual:

https://www.plumed.org/doc-master/user-doc/html/actionlist/

for some tips. Alternatively, if you want to have a chat drop me an email. I can't do this week (as I am at a conference) but I could probably find some time next week.

Good luck!

[PabloPiaggi]

Thanks, @gtribello !