How to define custom force in quantum chemistry packages?

Question

I don't know if this is really applicable but I want to add a custom force (say, $F(x,y,z)$) to some atoms in my geometry optimization. Is there any way to add such forces to the DFT or molecular mechanics calculations? (preferably GAMESS or Gaussian for QC or Avogadro for MM).

It can be also done as Molecular Mechanics (MM) as I just need an approximate geometry after applying the force.

Let me give an example. Say I have molecule A and molecule B and a force that changes in space as $F(x,y,z)$ but is only applied to molecule B. Now I want to calculate the new geometry of A and B.

I can think this is like adding a $V_F$ to the external energy in the Schrödinger equation, but I don't know how to do that in QC packages:

$$\left[-\frac{1}{2}∇^2 + V(r) +\color{red}{V_F(r)}\right]ψ(r) = Eψ(r)$$

Update: The idea of making a unit cell and making shrinkage on it looks great (suggested by Geoff Hutchinson here), but the problem is I couldn't figure out how to install CP2K on MacOS and also I am looking for something with a simple interface like Avogadro. I even tried Maestro by Schrödinger software, but couldn't make unit cell shrinkage for that.

Answer 1

The short answer is not. You cannot define a constant force for an atom or set of atoms in a system.

If you still want to modify force constants, this is the way to do it:

1. Gaussian, you can modify through the keyword GEOM(IHarmonic=n), where n is part of: $\frac{n}{10^6}$ which indicates the respective force applied to the system. Its units are $\frac{Hartree}{Bohr^2}$. To more info, please refer to this link.
2. If you want to introduce a specific force to an atom or set of atoms in MM, the best way to do it is by parameterizing your own system into your desire environment. The most important part is that you have to have in mind that the right part of the principal ec. of any MM describes the parameters of the specific force field of your desire system or set of systems under specific conditions (i.e. aminoacids in water like OPLS)

$$\vec{F}=\frac{\partial V}{\partial x}$$

3. Moreover, when you talk about forces in QM you are talking about the forces acting on electrons. Not on atoms.
4. But, if you still want to know a way to know how the forces of electrons (From DFT) and, at the same time, nuclei are, I recommend you Car-Parrinello. Please, refer to this link
5. If you have a particularly interest on simulating systems in QM with specific forces acting on electrons, the best way to do it (As far as I know) is to develop your own density functional.

Update

I Gaussian input file was requested to figure out how the keyword and input is used: Geom(iharmonic=n). I did not find the respective file, instead, I give a try twice by crating two input files with different values of n. The used Gaussian input files were:

chk=Test.chk
#p opt b3lyp/6-31g geom(iharmonic=n)

Test

0 1

H   
H 1 0.9 
O 2 1.2 1 120.

and the used values of n were 2 (On the left) and 20 (On the right). The respective Energy and RMS profile of both simulations are:

Also, in Gamess, you can use these keywords: Statpt, Force and its related keywords like Zmt. More info, Please see this link at page 2-121 and 2-131.

Answer 2

This is certainly possible in most MM packages and hence it's possible via QM/MM as well. I had to do a bit of digging to see how you might do this in each program (hint, it's not easy).

For Gaussian, your best bet is to read the extremely long section on molecular mechanics paramters.

For GAMESS, you'd need to install the Tinker program and use the QM/MM integration.

Avogadro and Open Babel don't currently allow custom forces, although it wouldn’t be too hard to code.

I think the other alternative, depending on why you need the forces, is to use a geometry constraint during optimization. For example, ensuring certain distances are retained.

The problem is basically, that without care, your geometry optimization will be forced to use numeric derivatives, since the QM or MM program won't know much about the functional form of your custom force.

Answer 3

If you are wanting the capability to apply forces to nuclei, to alter the equilibrium geometry, I do not know of any software with "a simple user interface like Avogadro" that offers built-in functionality to do what you seek. You are probably facing either a programming/scripting project or a great deal of manual iteration.

My recommendation would be that you investigate ORCA. ORCA incorporates the ability to perform geometry optimizations on systems where energy and gradient information is computed externally (ExtOpt mode; see Section 6.13.5 of the ORCA v3.0.3 Manual, downloadable from the forum linked above). I am actually planning to develop something rather like this, though with functionality not nearly as general as what you describe [arbitrary $F(x,y,z)$]. I plan to write a Python application that calls ORCA in this ExtOpt mode, and then within each ExtOpt cycle:

1. Makes a normal call to ORCA to compute the single-point energy and gradient in the absence of the external force field
2. Calculates the contributions of the external force field to the energy and gradient
3. Merges the QM energy and gradient computed by ORCA with the force field energy and gradient
4. Passes the composite energy/gradient back to ExtOpt ORCA

ORCA supports semi-empirical methods as well as full quantum methods, so even large systems should be tractable in a reasonable time frame, once you have built the execution framework to carry it out.

EDIT: For your application, where you're trying to simulate increased pressure, the 'constrained volume' approach would probably work pretty well using a piecewise potential function, where $V = 0$ within the molecular volume calculated for a particular pressure, and then ramps to a very high value at the borders of the box. The task then becomes appropriate definition of the box boundary, a problem which I suspect is similar to that involved in defining the cavity in implicit solvation models (COSMO$^1$, for example).

$^1$ Klamt and Schüürmann. Perkins Transactions 2 2: 799-805 (1993). doi: 10.1039/P29930000799