I have a technical question about using EMC scripting to implement a custom force field on a carved crystalline structure. My goal is to construct a spherical SiO2 silica nanoparticle using known force field parameters and “smooth” out its surface to consist of oxygens and silicon atoms that are coordinated with three oxygen atoms (instead of four), which I will then functionalize.
So far, I have an emc script that defines a rectangular unit cell, alpha_quartz_r.emc, and a script that imports it and builds a large crystal then carves out a spherical nanoparticle of a user-defined radius, build.emc. While my unit cell file includes the force field parameters that I need, I have also created a pair of emc-style force field files that contain the atom types, parameters, and rules that I will need, silica.prm and silica.top.
My question is as follows: is it feasible to, after performing my build step, import the field I have defined into my build script, and then reapply it to change the types of my existing atoms? I envision that I can “clean up” the surface in the way I described above by, for example, making all silicon atoms that are only connected to one or two oxygens a particular atom-type or group and then deleting all members of that group. This will also help me fix the issue of there not being charge balance in the system after carving. When I try to use the field = {mode -> apply} command after my script, I end up with the following output:
Warning: skipping group name 'obulk'.
Warning: skipping group name 'sibulk'.
Warning: no rules for group 'OBULK' in any field.
Warning: no rules for group 'SIBULK' in any field.
Error: core/fields.c:433 FieldsApply:
Missing rules.
Program aborted.
Even though, I do have rules for the SMILES that are associated with those group in the field. This makes me think that somehow, I have to reference the external field file in my unit cell file, but I’m not quite sure how that’d be done. Do you know if there is a better way to define your types and groups in the .emc file if you know you will be importing a force field file for atom typing?
I applaud you trying to use EMC Script. I tend to avoid it, when I can. We normally use EMC Setup for building around import structures. For your problem, such a script would look as follows:
# Options section
ITEM OPTIONS
replace true
field charmm/iff/clay, charmm/c36a/water_ions
number true
density 1
direction y
build_dir .
shape 2
emc_execute true
ITEM END # OPTIONS
# Groups section
ITEM GROUPS
water:f=water_ions O
ITEM END # GROUPS
# Clusters section
ITEM CLUSTERS
molecule import, name="SiO2_4H", &
type=surface, mode=insight, field=clay, percolate=true, &
tighten=1.5
water water, 1000
ITEM END # CLUSTERS
Here, I am using the IFF clay force field, which is an averaged version of the original version available online. I chatted with Prof. Heinz and he agreed with this approach. This particular structure has a y-direction orientation, hence direction y is set under options. I use the InsightII format for the unit cell, since it allows for defining crystal connections. The imported structure is a surface, as defined by type. I want to specifically use the IFF force field for this structure, as is specified by field. I need to set percolate=true since the structure has connectivity through its periodic boundaries. EMC Setup will automatically determine the size of the supercell to accommodate shape 2 under options. Supercells can be changed manually by setting ncells during the import. The option tighten=1.5 tightens the imported structure down to 1.5 Angstrom into the growth direction.
Additionally I am including 1000 water molecules to be typed with the CHARMM water_ions force field.
I am attaching setup.esh, the silicon oxide unit cell (SiO2_4H.car and SiO2_4H.mdf), and the IFF clay force field (clay.prm and clay.top). Place the clay force field in the directory $EMC_ROOT/field/charmm/iff/ in order for the above setup script to work. That way you can also use it in different projects.
Thank you so much for the response - I was thinking that likely importing the pre-built structure into a new script with the custom field would be the right way to go. I will take a closer look at this and get back to you.
As for using EMC Script… I will admit it is a bit confusing, though it seems there is some functionality that is very powerful, especially things like the carve command. Do you know if, for for example, that command could be accessed from an EMC Setup file? I have seen a couple of examples on here that use ITEM EMC to include scripting commands in .esh files.
I am getting closer to figuring this out. I have tried importing by built structure, aq, in three different ways: as a .emc file, as a .pdf/.psf, and as a .mdf/.car file. I am attaching all five of these files here. Some issues I am encountering:
→ attempting to retype the atoms when importing the .emc file leads to only the groups originally defined in the build script being able to be retyped, hence the issue is the same because all O or Si atoms end up being assigned to the same original type.
→ when I try to export the InsightII files for my spherical structure from EMC, the resulting .mdf file does not include connectivity information (only the atom types of the connected atoms), so it can’t be imported to EMC Setup.
→ perhaps the closest I have gotten: when I try to import the .pdb file, it actually works - but the field typing fails because, for some reason, erroneous formal charges are being assigned to the SMILES strings that are derived from the structure. I have tried using the formal=false flag in my import, but this does not change the behavior.
I am attaching the relevant files if you want to try to reproduce - I think that the pdb import will probably be the best bet if I can figure out how to fix the issue with formal charges. If all else fails, I may just try to write some kind of python script that deduces atom types based on the connectivity and retypes them accordingly.
Hi Pieter,
After further development- I have added some additional typing rules to get around the formal charge issue, and the atom typing after loading the .pdb/.psf files seems to work well for the bulk. For the surface, I need to delete some atoms of a certain type that are not correctly coordinated, then reload and retype, until I have a cleaned-up surface.
I think I should be able to accomplish this in EMC setup using delete, such as with delete sites=si_g+o_g
in my ITEM OPTIONS section. However, I believe I have encountered a bug, as when I try to do so I end up with the following error: Undefined subroutine &main::oper called at /home/slayding/emc/v9.4.4/scripts/emc_setup.pl line 3148.
I have tried using various other keywords with delete, and get the same error every time, leading me to believe that this is a legitimate error.