Dear Karsu
If you want to replace atoms with chemical groups in an existing
structure, it's relatively simple to delete the atoms you don't want,
and paste coordinates for the new chemical groups or subunits into the
correct positions with the correct orientation. It's also easy to
connect them with bonds. I demonstrate how to do these steps (using
moltemplate) below (see steps 4-7).
However the more difficult issue is that you will be forced to
create 3-body and 4-body angle and dihedral (and improper)
interactions which join the old and new bonds (so the chemical groups
don't flop about in a crazy fashion). You have to do that whenever
you change the bond topology of the system. You can attempt to
manually add these interactions to the "Angles" and "Dihedrals"
sections of the data file yourself, but there are often a large number
of them, and it is easy to make a mistake and omit one.
---- Figure out the new angles and dihderals manually ---
If you have figured out all of the old old and new angle and
dihedral interactions (to replace a C=O bond with a C-CH3 bond, for
example), you can define molecule templates before and after the
reaction and use "fix bond/react" to perform the reaction immediately
(setting Nevery to 1 and the reaction probability to 1).
https://lammps.sandia.gov/doc/fix_bond_react.html
Fix bond/react is powerful, but I'm not aware that it can change the
number of atoms in the system yet. (Hopefully somebody will correct
me if I'm wrong. This is relevant since you will be replacing O atoms
with CH3 groups.) Either way, I suppose you could delete the terminal
oxygens and place the methyl groups in the appropriate locations and
wait for the reaction to happen (which should hopefully take only 1
Verlet iteration).
---- Use a molecule builder to generate the new angles, dihedrals,
impropers ----
Alternatively, molecule builders allow you to specify the list of
bonds, angles, dihedrals that you want exactly, as well as
(optionally) a set of force-field rules which tell the software how
you want it to create these angular interactions automatically so you
don't have to.
If you are going to use a molecule builder, you have two choices:
A) Familiarize yourself with the file format the software uses to
describe these rules, so that you can tell the molecule builder how
you want it to generate new angle and dihedral interactions when you
replace an O atom with a methyl group (for example).
B) Choose an existing popular force-field (and choose a
molecule-builder which supports it), and change all of the atom types
in your structure to be compatible with the atom types used by that
force field.
Details (B):
Choose an existing popular force-field (and choose a
molecule-builder which supports it), and change all of the atom types
in your structure to be compatible with the atom types used by that
force field. Most molecule builders use atom type (and sometimes bond
type) to lookup angle, dihedral, and improper interactions, so you
will have to change the name of these atom types to something that the
force field understands. (You will also have to delete all of the
existing text in the "Angles", "Dihedrals", and "Impropers" section of
your DATA file, since all of these interactions will have to be
generated from scratch. It's not a great idea to combine interactions
from two different force fields in same simulation, or same molecule)
Then, add the new atoms (and if necessary, bonds). The molecule
builder will do the rest. There are instructions for doing this with
moltemplate below. (moltemplate currently supports OPLSAA, GAFF,
GAFF2, COMPASS-published.
---------- instructions for doing this with moltemplate ----------
---- step 1 ----
Firstly, I assume that you already have a LAMMPS data file (and the
corresponding input file that works with it, which I named "FILE.data"
and "FILE.in" in the example below). If so, you can create a
moltemplate object describing it using the "ltemplify.py" file
converter (which is included with moltemplate). This will enable you
to edit and modify your original system with moltemplate:
ltemplify.py -name OldSystem FILE.in FILE.data > old_system.lt
(See appendix B of the moltemplate manual if you want more information
about what "ltemplify.py" does.)
---- step 2(A) ----
Note: There are two versions of this step.
Either choose step 2A OR step 2B, but not both.
IF you decide to create rules to add new angle and dihedral
interactions to your system, then create these rules in a file (eg
"force_field.lt")
ForceField {
#rules and force field parameters go here...
}
Details: The file format of "force_field.lt" is described in the
"force_field.lt" and "oplsaa_simple.lt" files included (with context)
in these examples:
http://moltemplate.org/examples/force_fields/
A gruesomely detailed explanation of the file format is included
in slides #51-54 of this talk:
http://moltemplate.org/doc/moltemplate_talk_2013-8-07.pdf
...and it is also mentioned in chapter 6 of the moltemplate manual
(http://moltemplate.org/doc), as well If you have any questions about
the file format, let me know. (Incidentally, moltemplate will not
automatically erase your existing angle, dihedral, and improper
interactions. It will simply add new ones.)
Then insert the following line (import "force_field.lt") to the
BEGINNING of your "old_system.lt" file which you created in step 1, so
that the beginning of that file now looks like this:
import "force_field.lt"
OldSystem inherits ForceField {
---- step 2(B) ----
ALTERNATIVELY, if you don't want to create new force-field rules,
but you want to use the existing force-field rules from an existing
force field, you will need to replace all of the existing atom types
in your structure, with atom types from that force field (for example
OPLSAA). In that case, the procedure below is a little bit involved:
First, manually edit the "old_system.lt" file, and replace the first
line ("OldSystem {") with the following text:
import "oplsaa.lt" #(defines OPLSAA)
OldSystem inherits OPLSAA {
You will have to familiarize yourself with these atom types. For a
list of the atom type names used by the OPLSAA force-field that comes
with moltemplate, see this file:
https://github.com/jewettaij/moltemplate/blob/master/moltemplate/force_fields/oplsaa.lt
In particular, look at the write_once("In Charges") section, near the top.
Excerpt from that file:
set type @atom:80 charge -0.18 # "Alkane CH3-"
set type @atom:81 charge -0.12 # "Alkane -CH2-"
set type @atom:82 charge -0.06 # "Alkane >CH-"
set type @atom:83 charge -0.24 # "Methane CH4"
set type @atom:84 charge 0.0 # "Alkane >C<"
set type @atom:85 charge 0.06 # "Alkane H-C"
(Again, if you're curious about the file format simplified version of
this file is available here:
https://moltemplate.org/examples/force_fields/alkane_chain_single/oplsaa_simple.lt)
(Incidentally, if you're using OPLSAA, then avoid the atom-types from
1-56 whose descriptions end in "(UA)". Those are for coarse-grained
("united-atom") simulations. Don't use them in an all-atom sim.)
Now comes the difficult part: Read the file created in step 1
("old_system.lt"). In moltemplate, atom types are given names
beginning with "@atom:". The "ltemplify.py" script currently usually
generates atom type names such as "@atom:type1", but they might have
different names, depending on the contents of the DATA file you fed to
it in step 1.) You need to replace each of these old atom type names
(eg. "@atom:type6") with the corresponding OPLSAA-compatible name (eg
"@atom:80"). (Incidentally, "@atom:80" in OPLSAA corresponds the
the carbon atom from a methyl group.) One way to do that is to open a
text editor and manually replace all of these text strings yourself
("@atom:type6" -> "@atom80"). However that's probably a bit
cumbersome and error-prone.
Alternatively, you can copy the "oplsaa.lt" file into the working
directory containing your other .LT files. (such as "old_system.lt".
Incidentally, moltemplate will preferentially use the local version of
any .LT file.) Then edit this file and insert additional "replace"
commands for each of your atom types. Scroll down to the beginning of
the list of "replace" commands (it's a long file), and add commands
like this to the beginning of this list:
replace{ @atom:type1 @atom:83_b013_a013_d013_i013 }
replace{ @atom:type2 @atom:81_b013_a013_d013_i013 }
replace{ @atom:type3 @atom:84_b013_a013_d013_i013 }
replace{ @atom:type4 @atom:82_b013_a013_d013_i013 }
replace{ @atom:type5 @atom:85_b046_a046_d046_i046 }
replace{ @atom:type6 @atom:80_b013_a013_d013_i013 }
replace{ @atom:type7 @atom:78_b044_a044_d044_i044 }
replace{ @atom:type8 @atom:79_b045_a045_d045_i045 }
Add one such command for each of the atom types in your original data file.
What do these weird "80_b013_a013_d013_i01" strings mean?
Actually, I lied earlier. @atom:80" is not the full name of the
methyl carbon atom type in OPLSAA. (The full name is actually
"@atom:80_b013_a013_d013_i01", but nobody wants to type that in, so I
allow them to use the abbreviated name "@atom:80" instead. That's
actually the purpose of all of these "replace" commands.)
To recap, for each TYPE of atom in your "old_system.lt" file (eg
"@atom:type6"), you would need to:
FIRST, figure out the abbreviated atom type name in the force field
you want to use (eg "@atom:83")
THEN, figure out what the corresponding full-length version of that
atom type name (eg "@atom:80_b013_a013_d013_i01"), by searching for
@atom:80 in the list of "replace" commands in the force field file
("oplsaa.lt").
This is tricky, so you have any questions about this then email me.
Once this is done, then things get much easier:
---- step 3 -----
Create a new file (eg "system.lt") with the following contents:
---- system.lt ----
import "old_system.lt" #(defines "OldSystem")
old_system = new OldSystem #(create a single copy of the old system)
...Then, if you can write a script which can identify all of the
terminal oxygen atoms in your structure that you want to replace,
along with the atoms they are currently bonded to, then you can use a
list of moltemplate commands, added to the "system.lt" file, such as
these:
---- step 4 Delete the terminal oxygens ----
Note: The remaining steps assume you are using the OPLSAA force field.
If not, then replace OPLSAA with the name of whatever force field you
are using, and change the @atom type names accordingly.
OldSystem inherits OPLSAA {
# augment the definition of "OldSystem" with these commands
delete $atom:old_system/type7_243
delete $atom:old_system/type7_247
delete $atom:old_system/type7_253
}
:
# In moltemplate, anything beginning with $atom: (including the "/")
# refers to an "instance" of an atom instance (ie. an atomID counter).
# The name of the atoms may vary, but $atom names generated by
# ltemplify.py usually begin with "type?_" followed by a counter
# indicating which copy of that atom type you are interested in.
---- step 5: change atom types bound to oxygens ----
# Optional? Change the types of the atoms which
# were bound to the terminal oxygens:
# (eg "type7_151", 153, 156,...)
OldSystem inherits OPLSAA {
# augment the definition of "OldSystem" with these commands:
# overwrite the atoms which were bonded to the terminal oxygens:
write("Data Atoms") {
$atom:type7_151 $mol @atom:84 0.11 -3.365 0.678 1.133
$atom:type7_153 $mol @atom:84 0.11 -3.892 1.592 -0.905
$atom:type7_156 $mol @atom:84 0.11 -4.264 -1.466 -0.292
:
}
}
# Note: The coordinates must also be copied from the original data file.
# (unfortunately)
ch3group.lt (2.19 KB)