# How to define custom force in quantum chemistry packages?

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.

• Can you give a bit of insight into why you need the custom force - what you're trying to do here? – Geoff Hutchison Apr 19 '15 at 23:45
• Why are you trying to install CP2K on a laptop? Why don't you use a Linux cluster? – Geoff Hutchison Apr 20 '15 at 19:07
• @GeoffHutchison I tried it but as I said I am not familiar with Linux command line. It needs SVN which is already a big issue to get it work and then its installation is another problem itself. On the other hand I can not see any extension in Avogadro that can make an input file for CP2K. – Aug Apr 20 '15 at 19:25
• Are you wanting to apply forces to the electron cloud, or just to the nuclei? Or, both? – hBy2Py Apr 21 '15 at 14:08

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}$$

1. Moreover, when you talk about forces in QM you are talking about the forces acting on electrons. Not on atoms.
2. 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
3. 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.

• I may want to give a try to IHarmonic concept.hartree/bohr^2 dimension is energy/dimension^2 = force/dimension. It is not force and that's pressure. If you look at my comments with Geoff, pressure is exactly what I want to apply to the system. Can you give a little more explanation or even better a worked example of such Gaussian input file please? – Aug Apr 21 '15 at 15:19
• @Aug, thank you for commenting my post, please see the update of my post and have a good day :) – Another.Chemist Apr 21 '15 at 22:09
• Wow ! That's really amazing ! I have just a little confusion about the physical meaning of this initial harmonic. Is this equivalent to the external pressure ( as this external pressure is what I want to deal with ) ? If so, as I see 1 hartree/bohr^2= 1557 Pa. By setting n=20, we are dealing with 20 x1556x 10^-6= 0.031 Pa which has not really a physical justified meaning ( 1 atmosphere is 101325 Pa ) . It means even at 1 atm pressure we need to set n= 6.4 x 10^7 !! which is not really applicable . Is there anything I am doing wrong here ? – Aug Apr 21 '15 at 23:06
• First, Pascal is not a quantum unit. Second, you can find a description of the relationship between pressure, force and harmonics (From the quantum point of view) in here: physics.stackexchange.com/questions/147128/… See at the largest answer – Another.Chemist Apr 21 '15 at 23:19
• I thought this initial harmonic subject is worth to have a separate discussion. Would you have a look at it please : chemistry.stackexchange.com/questions/29143/… – Aug Apr 22 '15 at 2:31

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.

• Thank you Geoff. Actually I want to simulate the behavior of some gas compound materials under pressure. Here I am trying to add the total energy of P.ΔV to the system. First I thought I may use the idea of freezing the distance of some atoms as you said. I should say the results are really promising and compatible with the experimental data but I am thinking when I want to publish the results, I may face critics as the whole system may not collapse uniformly under the pressure. Thats why I want to do something more sophisticated. – Aug Apr 20 '15 at 0:05
• If you want to simulate pressure, your best bet would be to use a code with periodic boundary conditions and change the unit cell size (i.e., volume). Since n, R, and T are constant, P and V are inversely related. – Geoff Hutchison Apr 20 '15 at 0:07
• Moreover, for such simulations, you might want to consider MD, e.g., using CP2K, which will let you perform both unit cell calculations and QM dynamics. – Geoff Hutchison Apr 20 '15 at 0:08
• I really like the idea of unit cell. Is it possible in Avogadro ?( I am more comfortable with Avogadro and also I have Mac and seems CP2K is just for linux ) . I tried making a larger unit cell from crystallography menu, ran UFF for two methane molecules as a test and then tried to make the unit cell smaller . but the whole molecular system burst ! – Aug Apr 20 '15 at 0:31
• At the moment, the UFF optimizations in Avogadro and Open Babel do not respect unit cells. You can, however, use Avogadro to create the unit cell and then use CP2K or other program for the simulation. For publication-quality results, that would be my suggestion (i.e., LAMMPS, Tinker, GROMACS, CP2K, etc.) – Geoff Hutchison Apr 20 '15 at 0:52

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