Force field (Dreading, UFF) do not work when space group is used for structure descrition

GULP
1

I had created GULP input files by Materials Studio (Materials Studio stated GULP supports calculation only in P1 space group). The generated .gis works OK with Materials Studio independent Windows GULP build …

If I modify the files to describe the structure by the help of space group (tested P21, P-1) the geometry optimization generate nonsense …

Any idea what is wrong ?
The nonsence generating .gis follows (DMSO structure) :

opti conv conse qok nomod pres noauto
ftol 1e-005
gtol 0.0001
xtol 1e-005
maxcyc 1000

title
GULP calculation from Materials Studio for DMSO
end

cell
5.173000 5.636000 6.883000 71.510000 83.800000 63.070000 1 1 1 1 1 1

fractional
S1 core 0.601900 0.164000 0.149900 0.000000 1.000000 0.0 1 1 1
O1 core 0.330000 0.121000 0.203600 0.000000 1.000000 0.0 1 1 1
C1 core 0.907000 -0.207000 0.248000 0.000000 1.000000 0.0 1 1 1
H1 core 0.914000 -0.313300 0.161000 0.000000 1.000000 0.0 1 1 1
H1 core 0.872300 -0.293900 0.385400 0.000000 1.000000 0.0 1 1 1
H1 core 1.089400 -0.199100 0.245500 0.000000 1.000000 0.0 1 1 1
C1 core 0.666000 0.284000 0.337000 0.000000 1.000000 0.0 1 1 1
H1 core 0.528300 0.475900 0.315800 0.000000 1.000000 0.0 1 1 1
H1 core 0.859000 0.268900 0.326400 0.000000 1.000000 0.0 1 1 1
H1 core 0.645200 0.173500 0.470400 0.000000 1.000000 0.0 1 1 1

Species
C1 core C_3
H1 core H_
O1 core O_2
S1 core S_3

spacegroup
P -1

connect 1 2 double 0 0 0
connect 1 3 single 0 0 0
connect 1 7 single 0 0 0
connect 3 4 single 0 0 0
connect 3 5 single 0 0 0
connect 3 6 single 0 0 0
connect 7 8 single 0 0 0
connect 7 9 single 0 0 0
connect 7 10 single 0 0 0

library dreiding

output cif DMSO-out

2

Hello Michal,
Obviously, the connectivity is not set correctly. Add the “mole” keyword, remove “noauto” from the input line, and also remove all the “connect” directives. They seem to be incomplete. I checked your input using gdis and certain bonds appear broken. The input I suggest above, works fine for me (Gulp 6.2 compiled on windows11 by mingwin/msys). Also, it might be better to add charges.

3

The input was working OK in P1. I had only removed the symmetry generated second molecule and changed space group from P1 to P-1. So it is impossible the input is wrong.
The connectivity is correct for 100%, it can be checked easy manualy based on the DMSO formula.
How can I remove “connect” directives from the file, when they gives the code the instruction witch forcefield constants should be used ?
GULP is not able to identify atom types and bond types for forcecfields automaticaly - this must be specified.
Can you, please, post the input witch works for you ?

4

The following works fine. The atom types are already given in the species section, and the bonds are inferred by gulp using the “molecule” keyword.

opti conv mole conse qok nomod pres
ftol 1e-005
gtol 0.0001
xtol 1e-005
maxcyc 1000

title
GULP calculation from Materials Studio for DMSO
end

cell
5.173000 5.636000 6.883000 71.510000 83.800000 63.070000 1 1 1 1 1 1

fractional
S1 core 0.601900 0.164000 0.149900 0.000000 1.000000 0.0 1 1 1
O1 core 0.330000 0.121000 0.203600 0.000000 1.000000 0.0 1 1 1
C1 core 0.907000 -0.207000 0.248000 0.000000 1.000000 0.0 1 1 1
H1 core 0.914000 -0.313300 0.161000 0.000000 1.000000 0.0 1 1 1
H1 core 0.872300 -0.293900 0.385400 0.000000 1.000000 0.0 1 1 1
H1 core 1.089400 -0.199100 0.245500 0.000000 1.000000 0.0 1 1 1
C1 core 0.666000 0.284000 0.337000 0.000000 1.000000 0.0 1 1 1
H1 core 0.528300 0.475900 0.315800 0.000000 1.000000 0.0 1 1 1
H1 core 0.859000 0.268900 0.326400 0.000000 1.000000 0.0 1 1 1
H1 core 0.645200 0.173500 0.470400 0.000000 1.000000 0.0 1 1 1

Species
C1 core C_3
H1 core H_
O1 core O_2
S1 core S_3

spacegroup
P -1

library C:"path to dreiding library file"\dreiding.lib (this modification may not be necessary in your case)

output cif DMSO-out

The above input produces reasonable results. Nevertheless, the use of P1 together with a defined connectivity may be necessary for the dreiding force field, due to the way it calculates torsions. I believe prof. Gale can clarify this.

5

image
image900×794 64.8 KB

This is what your original input produces in gdis. This is corrected by removing the “noauto” keyword. The “mole” keyword isn’t really necessary if you remove the “noauto”, and you can also keep the “connect” directives.

6

Thank you for the image. It sounds lik gdis or GULP do not work correctky with symetry operations. Insead of only applicating symmetry operations they in addition shift the symmetry generated atoms back to the unit cell. This should not hapen, becouse it breaks the connectivity. Such handlink is OK only for DFT calculations, not for force fields …
Thank you for the input file - I will check the results against reference Dreading implementation to see the force constats were correctly assigned.
I am almost sure your input creates S-O single bond instead of S=O so the result can not be correct (becouse automatic GULP bondings assign single bond to everithing).

7

I can confirm your input generates 100% identical results as reference Dreading calculation (done in Forcite+ code in Materials Studio, space group P1). It works in this case because Dreading force field does not identify the force-fields parameters based on bond type but based on atom types (so it ignores the single bond information and handle S_3-O_2 as S=O. I am not sure if this can work for more complex force fields when bond type (typically resonant) recognition will play a rule … But maybe resonant bond will be covered good enough by the atom type …
I still wonder what is GULP doing when it apply symmetry operation . For non covalent bonded inorganic structures it give sense to shift the symmetry generated atoms back to the 0-1 fractional coordinates range, but this is not the case of molecular crystals when bonds should be specified …

8

In that case, you can try the bondtype keyword to define the types of bonds precisely.