• Extract the centers of dimers from the MD trajectory
• Compute Q4/Q6 and check the separation
• In lammps, organize the molecular trajectories per dimer

https://pubs.acs.org/doi/10.1021/acs.molpharmaceut.1c00727

@Kalharap Here is the code to generate the cif file for the molecular centers

from pyxtal import pyxtal
for file in ['cifs/ACSALA.cif', 'cifs/ACSALA13.cif']:
    f1 = pyxtal(molecular=True)
    f1.from_seed(file, ['aspirin'])
    c1 = pyxtal()
    sites = [{}]
    s = f1.mol_sites[0]
    sites[0][str(s.wp.multiplicity)+s.wp.letter] = s.position
    c1.from_random(3, f1.group.number, ['C'], f1.numMols, lattice=f1.lattice, sites=sites)
    c1.to_file(file[:-4]+'-center.cif')

I tried to run the calculation to test the similarity CV with the following steps

• uploaded to cif files to obtain the environment pdb files from the online app
• followed the master class example to run two NPT MD simulation.
• computed the CV's histogram. The local folder contains all infos.

In the plumed.dat, I defined the CV as follows

ENVIRONMENTSIMILARITY ...
 SPECIES=1-1024:21
 SIGMA=replace
 CRYSTAL_STRUCTURE=CUSTOM
 REFERENCE_1=../env1.pdb 
 REFERENCE_2=../env2.pdb
 LABEL=s
 MEAN
... ENVIRONMENTSIMILARITY

• env1.pdb contains four environments from aspirin-I.
• env2.pdb contains four environments from aspirin-II.

However, I could not get a good separation.

(plumed-masterclass-2022) qzhu@cms GAFF (main) $ ./run.sh 
# Sigma  Overlap
0.01 0.9924291975
0.02 0.9482459949999998
0.03 0.8825642975
0.04 0.83610403
0.05 0.781628725
0.06 0.730893365
0.07 0.7065540025
0.08 0.7075560825
0.09 0.7220368050000001
0.10 0.7458443849999999
0.11 0.780875315
0.12 0.8296926125