Compute fep with atom id instead of types possible?

Dear all,

I am trying to calculate the free energy of adsorption of a molecule in a zeolite framework. I managed to adapt the examples given for the hydration free energy of methane to calculate the free energy of adsorption of methane in the zeolite.
Now I want to calculate the free energy of adsorption for different, more complex molecules, for example aspirin.
According to the documentation I need to perturb the charges of the molecule via
"compute FEP all fep ${TK} pair lj/cut/coul/long/soft lambda 12 34 v_dlambda atom charge 4 v_dqC atom charge 3 v_dqH & "
since lambda doesn’t affect the kspace term.
Now comes the problem:
As can be seen in the data file below, the atom types are determined by their respective LJ parameters and I have different charges for the same atom type. Therefore I cannot adapt the charges as compute fep inteds me to do by their respective atom type. I would like to adapt the charges by the ID of the respctive atom, but this does not seem possible.

Am I missing something here? Could someone maybe suggest a solution? Or does LAMMPS just expect all atoms of the same type to have the same charge?

I use parameters from the Sage forcefield generated by the OpenFF Toolkit, the datafile of the molecule looks something like this:
Aspirin parameters from Sage with OpenFF

21 atoms
21 bonds
32 angles
47 dihedrals
24 impropers

8 atom types
10 bond types
6 angle types
11 dihedral types
2 improper types
-3.8374 96.1626 xlo xhi
-3.5916 96.4084 ylo yhi
-1.4821 98.5179 zlo zhi
0.0 0.0 0.0 xy xz yz

Masses

1 15.99943
2 15.99943
3 15.99943
4 12.01078
5 12.01078
6 1.007947
7 1.007947
8 1.007947

Pair Coeffs

1 0.16846514 3.0251065
2 0.20947353 2.99716
3 0.2102061 3.0398122
4 0.086879315 3.4806469
5 0.10884061 3.3795318
6 0.015611343 2.5725815
7 0.015779483 2.6445434
8 1.2326e-05 0.53453923

Atoms

1 1 1 -0.3812 1.2333 0.554 0.7792
2 1 2 -0.6051 -0.6952 -2.7148 -0.7502
3 1 3 -0.565 0.7958 -2.1843 0.8685
4 1 3 -0.518 1.7813 0.8105 -1.4821
5 1 4 0.1801 -0.0857 0.6088 0.4403
6 1 4 -0.1476 -0.7927 -0.5515 0.1244
7 1 4 -0.147 -0.7288 1.8464 0.4133
8 1 4 -0.057 -2.1426 -0.4741 -0.2184
9 1 4 -0.082 -2.0787 1.9238 0.0706
10 1 4 -0.154 -2.7855 0.7636 -0.2453
11 1 4 0.6657 -0.1409 -1.8536 0.1477
12 1 4 0.6651 2.1094 0.6715 -0.3113
13 1 5 -0.1471 3.5305 0.5996 0.1635
14 1 6 0.154 -0.1851 2.7545 0.6593
15 1 6 0.159 -2.7247 -1.3605 -0.4564
16 1 6 0.141 -2.5797 2.8872 0.0506
17 1 6 0.143 -3.8374 0.8238 -0.509
18 1 7 0.0827 3.729 1.4184 0.8593
19 1 7 0.0827 4.2045 0.6969 -0.6924
20 1 7 0.0827 3.7105 -0.3659 0.6426
21 1 8 0.448 -0.2555 -3.5916 -0.7337

Thanks in advance for any help!

How about you define a group that contains only the atoms where you want to change the charge and then use compute fep only on this group instead of the group “all”. For the group definition you can use individual atom IDs.

That is a good idea, thank you!

Unfortunately another problem has arised, a bit more broad than the first one. I am trying to write an input, that can do free energy calculations on all of my systems with hardcoding as little as possible. As already mentioned, I want to get the free energy of adsorption or solvation of more complex molecules.

For the script to adapt the charges in the equilibration before the fep I iterate over all atoms of the molecule with code like this:

variable tt loop ${molecule_atoms}
label qt
variable id equal ${framework_atoms}+${tt}
variable adapt_q equal v_fwlmbd*v_def_q${tt}
set atom ${id} charge ${adapt_q}
next tt
jump SELF qt

In a similar way I can compute the charges, that compute fep perturbs to. Now for the compute fep command itself, I would like to be able to use a loop inside the compute command to iterate over the molecule atoms and specify the charges for the pertubation, so I can use the input for systems with molecules with various numbers of atoms. Similar to the first loop I tried:

variable ttt loop ${molecule_atoms}
compute FEP PPCP fep ${TK} &
    pair lj/cut/coul/long/soft lambda 1*2 3* v_dlambda &
label nqt &
variable id equal ${framework_atoms}+${ttt} &
variable adapt_q equal v_fwlmbd*v_def_q${ttt} &
variable pert_q${id} equal ${adapt_q} &
atom charge ${id} ${adapt_q} &
next ttt &
jump SELF nqt &

which LAMMPS doesn’t like, echoing Illegal optional keyword in compute fep.

So, my question is if using a loop inside a compute is even possible at all.

The (less elegant) alternative would certainly be to write a code, e. g. python skript, that writes the lammps input for each system of interest.

Kind regards and as always,
Thanks for your help!

That question is answered by the fact, that you get an error.

Parsing in general and thus also looping in LAMMPS input is line based. Using a continuation character ‘&’ collapses lines. But looping within a line is not possible.

Or you generate your input on the fly using the LAMMPS python module.

That is also a good idea! Thanks a lot!

Regarding your previous suggestion I am not sure if I understand you correctly or if I made myself clear.
Are you suggesting that by using

compute FEP molecule fep ${TK} ... 

I can perturb the properties of single atoms by referencing their id instead of whole atom types? If so, could you give an example of the syntax, since the documentation does not show such an option.

Or do you mean that if I use the group command I can exclude some atoms of the same type but with other charges if they are of the same type but not in the same group?
If you mean this, I think your suggestion does not work for me, because after all I need to include the whole molecule in the pertubation. Excluding atoms of the same type but with different charges would mean, leaving charges of the molecule, that is to be perturbed, untouched, right? This would not be my goal, because these atoms would not be included in the fep.

Please read my suggestion and the documentation for the fix command more carefully and pay attention to the fix group ID purpose. You can do what you want to do.

Are you suggesting to use the group command on each atom of the molecule, then making a compute fep for each group (namely each atom individually) and then use compute/reduce to sum the “per-atom”-pertubations?

No.

But you are suggesting to use the compute fep on a group of atoms in which each atom type only has one charge, right? So I can use the compute fep normally and perturb the charges by atom type?

No. Your thinking is far too convoluted and considering things that are impossible within what LAMMPS can do. But matters are usually simpler.

Unless the documentation for an individual compute style says differently, all computes “operate” on a group of atoms which is specified when you issue the compute command. That is why you must specify it. From the compute command docs:

compute ID group-ID style args

ID = user-assigned name for the computation
group-ID = ID of the group of atoms to perform the computation on
style = one of a list of possible style names (see below)
args = arguments used by a particular style

So you have two “selectors”, first the atoms that are within the specified group, e.g. your molecule and then e.g. the atom types you select on the command line, and compute fep would only operate on atoms that match both.

Yes, this is my understanding of the compute. But to me the problem still remains. Maybe I can illustrate my train of thought on a simple example and you could point me to where it derails?

Let me consider 1,1,1-Trifluoromethane in water. Let’s say the atom typing would put both C atoms in the same type, because their LJ parameters would be the same, but their charges would differ. The group “molecule” that I want to apply the pertubation to therefore consists of 8 atoms, and 3 atom types.
The compute command that I would use is:

compute FEP molecule fep ${temp} &
pair ... &
atom charge 1 v_pert-q1 &
atom charge 2 v_pert-q2 &
atom charge 3 v_pert-q3

but still I arrive at the point, where the charge of the atoms is set by their type and both C atoms would be perturbed to the same charge.
It seems to me both “selectors” are true, because both atoms are part of the molecule I want to perturb, meeting the first selector, and also have the same atom type, meeting the second one.

From my understanding, your solution would work for me only if I was excluding atoms from the group, for example only put the CH3-fragment in a group and apply the compute only to that group. Then I would avoid the problem but also I could not do the alchemical transformation in one step.

Thank you for the extensive help!

Then you can just assign a different, new atom type to those “conflicting” atoms and copy the LJ parameters over. When reading a data file, you can reserve space for additional atom types.

This would be a solution. The solution I chose was to edit the compute in the .cpp file so the pertubation is applied to atoms by their ID (or tag).
Is does seem to work now, even though some validation is still to be done.
Thanks again for the insightful discussion and the advice!