[lammps-users] compute the forces acting on a "dummy atom" and transferring them to nearby atoms

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps users.
I’m finding myself in need to solve a problem concerning a force-field for bisferrocene using a “dummy atom metod” described in
J. Phys. Chem. B 2006, 110, 33, 16640–16645.

I don’t know if a solution can be implemented within lammps.
The system is made by bisferrocene molecules added through the molecule command.

Each of the ciclopentadiene rings has a central dummy atom that binds to the central iron.
the dummy is a geometrical point computed at each step from the positions of the carbon atoms.
thus all the forces computed on the central atom at each step are transferred to the C1 and C2 atoms of the ring instead of directly acting as acceleration of the dummy atom and the position of the dummy atom is computed geometrically at each timestep from the position of the ciclopentadiene carbons.
I know nothing of cpp but I can handle python scripting

To sum up my problem is as follows:
I want to

  1. create x molecules of ferrocene with all the bonds/angles/dihedrals, including bonds linking the dummy atoms to the iron atoms and to the ring
  2. take the forces on the dummy atoms and add them to the fforces on the c1 and c2 atoms of the ring. then zero the forces on the dummy atom.
    3)integrate to compute the positions and velocity for the following timestep
    4)recompute the position of the dummy atom based on the position of the ring atoms.
    iterate. Is there a way to do this with a combination of compute and fix or is it something that requires external scripting? in the second case where do I place my script in order to have the operations occur in the intended order?
    Thank you in advance
    Mattia Siviero
2

Mattia,

I don’t see a simple and clean way to implement this without significant C++ programming. The primary way to have such dummy atoms implemented in LAMMPS (and I assume there are only non-bonded interactions with it) would be to compute its position and project the forces on it back to the adjacent atoms as part of the pair style.

The only approach that doesn’t involve programming would be to make the cp-rings including the dummy atom (which must be given finite mass but it can be a very small one) rigid objects (everything else can be flexible).

Anything involving fixes and computes would result in positions and forces being “out of step”.

Axel.

3

Thanks for the quick replies
To be fair my attempt at modeling ferrocene is part of a larger problem which is trying to drop coat a gold substrate with ferrocene terminated thiols. I tried to first work without centroid, but any such approach risks to overconstrain the problem. In practice I had issues determining force constants that were just right and prevented the molecule from blowing up. So I was trying to find a way to make centroid work without having them escape the center of the ring.

I only have one doubt: making the whole can making the whole cp rings rigid cause the non rigid bonds and the translation of the center of mass to have higher energy? and most importantly, does this really affect the validity of the model?

4

I don’t quite understand what you are asking.

Considering the stability and properties of ferrorcene, making the cp rings rigid is quite a mild approximation. I would consider it similar to other common approximations like keeping bonds involving hydrogens constrained. But only rigorous testing can confirm it.

Axel.

5

Thank you for the clarification. I’m looking at the details of fix rigid. I would like to check my understanding of it:
I have to create the molecule file specifying very stiff bonds/angles/dihedrals in order for energy minimization to preserve the desired structure.
I have to add the subgroup of atom I want to become a rigid body to a group.
Finally I have to apply fix rigid with the “molecule” style. This will treat the atoms of each molecule that are in the fix group as an independent rigid body linked to the molecule by all existing bond/angle/dihedral styles that were in place before fix rigid was applied. and it will leave the other atoms free to move withing the limits of their bonds.
Is that correct?
thank you,

Mattia Siviero

6

Thank you for the clarification. I’m looking at the details of fix rigid. I would like to check my understanding of it:

depending on the specifics of your system and needs, you may want to use one of the other rigid fixes, possibly the “small” variant.

I have to create the molecule file specifying very stiff bonds/angles/dihedrals in order for energy minimization to preserve the desired structure.

since you don’t have complete force field settings for the entire ferrocene structure, it is probably best to avoid minimization altogether and use simulated annealing for equivalent purposes.

I have to add the subgroup of atom I want to become a rigid body to a group.

Finally I have to apply fix rigid with the “molecule” style. This will treat the atoms of each molecule that are in the fix group as an independent rigid body linked to the molecule by all existing bond/angle/dihedral styles that were in place before fix rigid was applied. and it will leave the other atoms free to move withing the limits of their bonds.
Is that correct?

this is all explained in the documentation and if you have doubts, you should make tests (I would do those anyway) with small subsystems to verify that everything is working as expected. that is good practice in any case. even if I agree with what you describe, that does not automatically mean that you are implementing everything correctly. there are always lots of ways to make mistakes and it is impossible to predict them from remote.

ciao,
axel.