Lattice Energy Optimization with pGFN-FF

Hello, I have been trying to use the new pGFNFF model to optimize for lattice energies of zeolites using GULP. I just started trying to use the gfnff keyword with the gwolf keyword along with opti and conp to compare to a run using the same crystal structure via the core-shell method that a colleague is using. We use the same opti and conp keywords, but the value I get for total lattice energy from GULP is approximately 1/10th the magnitude of what my colleague gets from his test using the core-shell method. I have tested vs 2 other zeolite structures and had similar results in my lattice energy maintaining a 1/10th magnitude. As I am new to using this program I am more than likely making a few mistakes in my input file or there is something I am missing in general. Attached are the input files I am using to compare along with terminal outputs for each run as well. Any help/clarity on what is going on would be appreciated.
Test_Run_ABW0.txt (26.5 KB)

ABW0.gin (1.1 KB)
ABW0.gout (3.7 KB)

Hello. This is not actually a case of a code issue, but an understanding of force fields problem. With force field models (unlike all-electron quantum mechanics) the absolute energy of a structure doesn’t mean anything since it depends on the choices you make. A core-shell model will often give large energies due to formal charges, whereas partial charge models will give lower total energies. This doesn’t make one better or worse - they’re just different. What matters are the energy differences, structures, properties, etc.

Hello Dr. Gale this was originally my question that I had to ask my friend Emmanuel to post for me because I lacked the ability to attach files to the forum because this is my first time posting. I am an undergraduate student trying to assist in research, and as you pointed out I very much do lack a solid understanding of force field models. If possible could I ask if you have any of the documentation for using GFN-FF in GULP or advice as to which keywords I ought to be using or not. Thank you for your patience and assistance, it is greatly appreciated!

Hi Bob. Thanks for the explanation. For keywords, all you should normally use to activate this force field is “gfnff gwolf”. The first is absolutely necessary, while the second is my personal recommendation. In terms of more details of GFN-FF and the periodic implementation then the place to go is the paper:

Gale, J.D., L. LeBlanc, P.R. Spackman, A. Silvestri and P. Raiteri, “A universal force field for materials, Periodic GFN-FF: Implementation and examination”(https://pubs.acs.org/doi/10.1021/acs.jctc.1c00832)
J. Chem. Theory Comput., 17, 7827-7849 (2021)

Plus of course the original paper of Spicher and Grimme in Angewandte Chemie.
Regards,
Julian

Sorry to bother you again Dr. Gale, but I have looked at your paper, and specifically looked at the supplemental information. I have been looking at table S2, and trying to reproduce that information using GULP and the cif file from the CSD database that I convert to a .gin file with the gfnff conp opti keywords to compare to the dATM row for both acetic acid and ammonia (I noticed that by default pgfnff seems to dampen the ATM 3-body dispersion term enabled). Am I correct in my understanding that, that table was filled using the current implementation of gfn-ff in GULP? When I try to run the calculation I end up with different answers for both the unit cell parameters and the lattice energy. Thank you.

I figured these files might be of help with being able to point me in the correct direction.

ammonia0.gin (750 Bytes)
ammonia0.gout (2.6 KB)
acetic_acid0.gin (1.4 KB)
acetic_acid0.gout (4.8 KB)

Thanks again.

Hi Bob,
When I run with your input I get a lattice parameter of 4.82 Ang, versus 4.83 Ang in the SI for the same settings. This is pretty close and the small difference between the two may be because you’re not using symmetry or other minor numerical differences. For the lattice energy, your input gives -71.4 kJ/mol (to 1 d.p.) which is exactly the same as the SI. My guess here is that you’ve misinterpreted the “lattice energy” that GULP prints out. This is the total energy for the system in the case of GFN-FF. To get to the standard lattice energy you have to divide by the number of molecules per unit cell (4) and subtract the intramolecular energy of NH3 (which can be obtained from a gas phase calculation). Putting this together and converting from eV to kJ/mol you should get the same answer as in the SI.

Note, I’d recommend not using TAB characters in the input as a separator, though they don’t seem to cause issues here.

Thank you for all your help, that was indeed the issue. I had divided by 4, but forgot about subtracting the energy of the molecule itself. Thank you for also pointing out the potential issues with using the tab character.