Rigid molecules with GULP 6.0

Dear Dr. Gale,

I am interested in optimizing structures that contain water molecules coordinated to the metal centers of a periodic framework. I started testing the usual models (tip3p, tip4p, and so on), that are rigid. Thus, my system will contain a flexible framework allowed to relax, as well as the position of the waters. However, waters’ bonded terms should be kept as defined by their force field.

I saw the new feature of ‘rigid molecule optimization with quarternions’ in GULP 6.0, but I could not set up correctly the input file. I’ve tried using the ‘molrigid’ definition, but I am not sure about how to properly define (Q1 Q2 Q3) values for the quarternions. Unfortunately, there is no example of an input using the ‘molrigid’ keyword. Besides, although I doubted it were good approaches, I have also tried setting different regions (cartesian/cartesian rigid; cartesian region 1/cartesian region 2 rigid); and defining harmonic potentials for the water bonded terms with really high k constants, unsuccessfully.

Considering my problem (optimization with flexible and rigid molecules):
a) what is the best/correct approach to declare this in the input? If it is the usage of ‘molrigid’ definitions, can you give me an example for properly defining it?

b) If we want to use a water model that contains dummy atoms (massless point charges), what is the correct definition of them? (I usually set an atom with an element different from the remaining ones, and change the values for the mass (to zero), as well as cova and vdw).

Any other comment about properly defining water molecules is really welcome.
Thanks in advance,
Best regards.

Hi Alexa,

The easiest approach to modelling your system would be to use a flexible water model, such as SPC/Fw for example, which is one that we use in a lot of our work. This avoids the complexity of using rigid molecules. However, if you really want to use rigid molecules then this would be possible. By far the easiest way to set things up is just by using the geometry of the molecule in the original structure (assuming you’ve adjusted the bond lengths and angles to be consistent with the force field). If you specify a molecule keyword and “rigid” then the code will take the geometry from any molecules that are found in the structure. This avoids having to use molrigid and define the quaternions yourself (this is more there for restarting of runs). Of course if you dump a restart file then these should be specified and so the geometry would have been converted into this format.
For dummy atoms (as required for TIP4P etc) then you are doing the right thing - just choose an element that is not in the system and then change the properties to be those that you want, such as zero mass).
Best regards,
Julian

1 Like

Thank you for your answer!
Best regards