converting com from unwrapped to scaled

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I’m working on a problem that involves finding the first hydration shell of a molecule I am interested in.

My approach to this involves finding the center of mass(COM) for molecule I am interested in and then using those coordinates as a new “molecule” and append the coordinates to a trajectory file to read into VMD to analyze further using the tools implemented in that program.

I have however ran into the problem that when I load the trajectory into VMD and display the new “molecule” or rather representation of the COM, I get a bead that is not really in the center of the molecules who COM I am interested in.

I have had this question before and still am unclear if I am addressing converting the center of mass into called units.

My python code is:

1 com1[j][0] = com1[j][0] - xl*m.floor(com1[j][0]/xl)
2 com1[j][0] = com1[j][0]/xl

which just takes the COM in unwrapped units and 1) determines how many box lengths to add or take away depending on the sign of the COM coordinate and 2) then divides by the box length to then get the scaled unit

ex. COM x -45 x box length: 51
after line 1: 6
after line 2: 0.117

I have checked my numbers I know the algorithm is working.

However I’m not sure if this is the correct algorithm as I don’t see physical results when I convert the COM and append them back onto the trajectory.

I print out the COM when I do the run. I also rerun the trajectory to get the COM by themselves to avoid a complicated code.
The two COM between the run and rerun match so I don’t believe the problem to be that.

Any thoughts on this problem is much appreciated

Thanks

Michael Goytia
Graduate Student
University of Utah Chemistry Department

2

Your line 1 makes a lot of assumptions. It might work for certain
atoms coordinates and simulation cell ranges, but it is not general.
Look at how LAMMPS handles this in domain.cpp:

    while (coord[0] < lo[0]) coord[0] += period[0];
    while (coord[0] >= hi[0]) coord[0] -= period[0];
    coord[0] = MAX(coord[0],lo[0]);

Aidan

3

Aidan,

Can you explain what this last part of the code

"coord[0] = MAX(coord[0],lo[0]);”

Thanks

Michael Goytia
Graduate Student
University of Utah Chemistry Department