Measuring pair/local distance of selected atoms

LAMMPS LAMMPS Beginners
1

I’m trying to measure and log the distance(s) anytime a calcium cation (2+) comes within the lj/charmmfsw/coul/long pair cutoff of two carboxylates (1- each) on side chains of aspartic and glutamic amino acids. Here’s an example of such an interaction:

image2
image2707×492 90.5 KB

this is my first try:

group oc type 39
group ca type 47
group ca_oc union ca oc

compute 1 ca_oc property/local ptype1 ptype2 patom1 patom2
compute 2 ca_oc pair/local dist
dump 2 ca_oc local 1 ca_oc.dump index c_1[1] c_1[2] c_1[3] c_1[4] c_2

since these interactions are fairly rare but the two carboxylate oxygens are always close to each other, and i want not to miss any, logging every 1 on a job running for 7 days creates a bloated dump file of several TBs which is ~98% O<->O.

(1) how can i tell dump local to only include interactions where ptype1 != ptype2 ?

(2) because charmmgui generates ions after tip3 water molecules, the atom IDs of the calcium ions are >1M which means dump local switches to E notation. Because dump local rounds to 5 decimals, i’m unable to distinguish which calcium ion is interacting with which carboxylate. what im actually trying to determine is which pair of carboxylates are interacting with a given calcium ion. is there a way to turn off exponential notation of dump local ?

Atoms

           1        1     33  -0.300    72.0175969762     3.6711606255   -13.6317739190 # PROA-1-MET-N-NH3
           2        1      7   0.330    71.2030649155     3.1557447367   -14.0297802498 # PROA-1-MET-HT1-HC
[...]
     1297527     2442     47   2.000   -32.3649286273   -59.3190680596    88.5585942907 # IONS-2442-CAL-CAL-CAL
     1297528     2443     47   2.000  -103.3356564158   -13.2219937221  -112.9722611123 # IONS-2443-CAL-CAL-CAL
     1297529     2444     47   2.000    61.4697979806    51.3687331027   -27.1765147868 # IONS-2444-CAL-CAL-CAL
     1297530     2445     47   2.000   -49.3672415527   -30.4155712039   -83.6149993609 # IONS-2445-CAL-CAL-CAL
     1297531     2446     47   2.000    80.6629122832   -40.9971418352   105.0700184218 # IONS-2446-CAL-CAL-CAL
     1297532     2447     47   2.000   105.8439127987   -41.0609458886    26.4601118012 # IONS-2447-CAL-CAL-CAL
     1297533     2448     47   2.000    -5.0039676348     6.9473668422   -21.5725091255 # IONS-2448-CAL-CAL-CAL
[...]

grep -v ’ 39 39 ’ ca_oc.dump:

ITEM: TIMESTEP
0
ITEM: NUMBER OF ENTRIES
528
ITEM: BOX BOUNDS pp pp pp
-1.2500000000000000e+02 1.2500000000000000e+02
-1.2500000000000000e+02 1.2500000000000000e+02
-1.2500000000000000e+02 1.2500000000000000e+02
ITEM: ENTRIES index c_1[1] c_1[2] c_1[3] c_1[4] c_2
26 47 39 1.29755e+06 10056 6.25891
27 47 39 1.29755e+06 10010 2.32465
28 47 39 1.29755e+06 10011 3.14261
29 47 39 1.29755e+06 10055 4.05748
30 47 39 1.29755e+06 10087 8.63207
31 47 39 1.29755e+06 10088 6.79561
32 47 39 1.29755e+06 9697 9.55493
33 47 39 1.29755e+06 9698 11.5767
266 47 39 1.29756e+06 3598 11.5345

(3) is there another better way to do this with a fix, a compute, or using the chunk functionality ? id prefer to measure the distance from a calcium ion to c.o.m. of carboxylate because its negative charge is delocalized between the oxygens due to resonance, not on the actual oxygen atoms.

2

Random thoughts:

Is this a job for compute rdf or compute adf?

You can “reorder” your atom IDs programmatically as follows:

  • Use delete_atoms and reset_atoms to write separate data files of protein, calcium, and water atoms
  • Use read_data with append to read in (for example) first protein, then calcium, and finally water data files to get the atom IDs in that order.

You can also use dump_modify format to change output formats of dump local columns.

It’s hard enough to do your desired processing “in-line” during the simulation; don’t even think about trying to track the residue COM. A suitable refinement is to track the carbon between each oxygen pair; that will be close enough to the centre of the delocalized double oxygen charges.

You can use compute coord/atom to track how many carboxylate O’s (or C’s) are near each calcium ion. You can use compute reduce ... max to determine if any calcium ion in the system has two or more neighbouring carboxylates, then pass that to dump custom with dump_modify skip to print out coordinates only on relevant steps. Print coordinates (as compared to pair interactions) because you can easily restrict printing by group (for example only carboxylate side chains and calcium ions) and limit yourself to N instead of about N^2 lines, whereas pair information with dump local is much harder to filter “in-line”.