Importing cyclodextrin (formerly nrings != nlinks error)

Dear all,

I plan to build a cyclodextrin model using EMC. However, a software error occurred during the modeling process, as follows:

Error: core/chemistry/ring.c:651 ChemistryRings:
       Found nrings != nlinks (7 != 8); increase scanning depth?
       Program aborted.

Even though I tried several other forms of SMILE structure, I still get the same error. My esh.file was written like this:

ITEM	OPTIONS

replace		true
field		pcff
density		0.01
ntotal		147

ITEM	END

# Groups

ITEM	GROUPS

mol		C([C@@H]1[C@@H]2[C@@H]([C@H]([C@H](O1)O[C@@H]3[C@H](O[C@@H]([C@@H]([C@H]3O)O)O[C@@H]4[C@H](O[C@@H]([C@@H]([C@H]4O)O)O[C@@H]5[C@H](O[C@@H]([C@@H]([C@H]5O)O)O[C@@H]6[C@H](O[C@@H]([C@@H]([C@H]6O)O)O[C@@H]7[C@H](O[C@@H]([C@@H]([C@H]7O)O)O[C@@H]8[C@H](O[C@H](O2)[C@@H]([C@H]8O)O)CO)CO)CO)CO)CO)CO)O)O)O

ITEM	END

# Clusters

ITEM	CLUSTERS

mol		mol,1

ITEM	END

How could I deal this problem, please?

Dear user,

Your chemistry is a carbohydrate super structure, which consists of seven six-membered rings and one overarching 42-membered ring, as shown in the following


EMC uses a scanning depth to recognize these rings, which by default is set to 8. In order to be able to build your system, you would need to increase this scanning depth to at least 42 – I would take 50 to be safe – which can be done by adding emc_depth 50 to your ITEM OPTIONS block. The resulting .esh script would read

#!/usr/bin/env emc_setup.pl

ITEM	OPTIONS

replace		true
field		pcff
density		0.01
number		true
emc_depth	50
emc_execute	true

ITEM	END

# Groups

ITEM	GROUPS

mol		C([C@@H]1[C@@H]2[C@@H]([C@H]([C@H](O1)O[C@@H]3[C@H](O[C@@H]([C@@H]([C@H]3O)O)O[C@@H]4[C@H](O[C@@H]([C@@H]([C@H]4O)O)O[C@@H]5[C@H](O[C@@H]([C@@H]([C@H]5O)O)O[C@@H]6[C@H](O[C@@H]([C@@H]([C@H]6O)O)O[C@@H]7[C@H](O[C@@H]([C@@H]([C@H]7O)O)O[C@@H]8[C@H](O[C@H](O2)[C@@H]([C@H]8O)O)CO)CO)CO)CO)CO)CO)O)O)O

ITEM	END

# Clusters

ITEM	CLUSTERS

mol		mol,1

ITEM	END

# EMC

ITEM	EMC	phase=1	spot=2

test_chirality	= {
  mode		-> true
};

ITEM	END

I added number true to avoid having to provide ntotal. This option means, that the number at the end of mol mol,1 is interpreted as the number of molecules or clusters to build. I also added a verbatim EMC Script paragraph to test chirality of your structure. This paragraph is added right after the build in build.emc and produces the stereo chemical state. It allows you to see, which stereo chemical centers stay the same, and which change after multiple builds.

A word of caution with respect to stereo chemistry: EMC produces correct stereo chemistry for linear molecules. However, it makes mistakes when considering ring structures like sugars when the ring closing atom is concerned. This ring closing atom is defined by your SMILES. I am aware of the issue and have a route to solution, but unfortunately haven’t had the time to fix it, since the fix is rather involved. Please decide yourself, if correct stereo chemistry is needed for the model you are developing. In the above case, I would roughly estimate, that 75% of the stereo chemistry is correct.

The above example can also end up not being build. This happens either because ends are not found or the molecule is growing too much into itself. Three settings influencing the build procedure are niterations, radius, and nrelax, of the which the defaults are 1000, 5, and 100 respectively. The niterations option sets the number of trials to find a solution before aborting. One could try to increase the number of iterations. The radius option controls the radius within atoms are relaxed. The center of the sphere is the added atom. I would not go higher that 10 Angstrom. nrelax controls the number of relaxation cycles EMC uses. Upping this number will result in more relaxed structures. Increasing either radius or nrelax or both will slow down the building.

Dear Prof. Veld

Thanks to your quick and detailed reply, the modified code was able to generate the model. However, as you said, the stereochemistry of the generated cyclodextrin is incorrect, and the chirality of the D-glucose that makes up the cyclodextrin is not consistent.

Do you think there is a way to solve this problem? Or I sincerely ask you, do you think there is any way to automatically assign a pcff force field to a manually plotted cyclodextrin if the atomic coordinates of the cyclodextrin are plotted manually first?

Dear user,

I am aware of the issue in EMC and should be fixing it. I know how, but unfortunately have not found the time to do so. A work around for the study of regular structures would be to build one structure in EMC, after which you manually correct the mistakes – in your case four in total. You can then import this structure by

molecule	import,name="molecule",type=structure,mode=pdb, &
			tighten=3,ncells=nx:ny:nz

where molecule represents the name of your imported structure and ncells controls the number of structure copies. The default is nx=1, ny=1, and nz=1. Say you would like to use eight structures. You could accomplish this by replacing nx, ny, and nz by the number 2. You would add a water group and cluster for embedding these structures in liquid water. Be aware, that you have to specify enough material to solubilize your cyclodextrin (as is the case for How to embed a carbon nanotube in a polyethylene melt).

Dear Prof. Veld

Thank you for your suggestion.

I have drawn the .pdb file of cyclodextrin like this (mol.pdb (22.0 KB)).
But I don’t have a .psf file. I would like to ask is it possible to generate .psf file in EMC based on this pdb file I have?

If I just import this .pdb file like the following,

ITEM	OPTIONS

replace		true
field		pcff
number		true
density		0.1
emc_execute	true

ITEM	END

# Clusters

ITEM	CLUSTERS

molecule	import,name="cd",type=structure,mode=pdb,tighten=3

ITEM	END

I will get the following error,

Error: core/format.c:359 FormatOpen:
       Error opening file "cd.psf".
       Program aborted.

Dear user,

The error you are experiencing tells you, that you are providing a .pdb, but not a .psf file. You would need to provide the file in order for EMC to be able to interpret your structure correctly. EMC needs a .psf for correct connectivity, since a .pdb does not contain this information (unless your using PDB the CONECT keyword, but for EMC this is a more experimental approach). I used VMD to produce a .psf. The resulting structure is given by

single.pdb (11.4 KB)
single.psf (14.1 KB)

which in combination with

ITEM	OPTIONS

replace		true
field		charmm/c36a/carb, charmm/c36a/cgenff, charmm/c36a/water_ions
ntotal		1000
density		1
emc_execute	true

ITEM	END

ITEM	GROUPS

water:f=water_ions	O

ITEM	END

ITEM	CLUSTERS

molecule	import, name="single", type=structure, mode=pdb, &
			tighten=3, field=carb, depth=64
water		water,1

ITEM	END

results in

Please note that you will have to download a new version of EMC (still v9.4.4 of 20230801, just a newer updated version, to be downloaded from here), since I noticed a bug existed in the C code.

Dear Prof. Veld,

I referred to your method for modeling. But I’m very sorry that I’m having trouble exporting the psf file using VMD.
May I ask which function of VMD you used to export psf file, as I was not able to export pdb file to psf file through VMD.

Dear user,

You can generate a PSF within VMD by means of starting vmd with vmd -dispdev text -f mol.pdb, after which you execute the following within VMD

set atoms [atomselect top all]
$atoms writepsf mol.psf
exit

Note, that VMD will warn you about the absence of angle, dihedrals, improper, and cross terms. However, the resulting mol.psf will include the necessary bond section, since your original mol.pdb included bond information through CONECT entries.

Dear Prof. Veld,

Thank you very much for your detailed explanations! With your suggestions, I can successfully build my target model now.