esp restrain

Whether it is possible by using the keyword «restrain» of «esp» group indicate that I want to draw the partial charges not to zero, but to other values I specified?

The options in constrains (just above restraints in the documentation) should be able to allow you to do what you want I believe.

Bert

The options in constrains (just above restraints in the documentation) should be able to allow you to do what you want I believe.

Bert

I have no clue what RESPA stands for and can't find anything online about it either.

I guess you want to tell ESP to give an atom a minimum or maximum charge. This is not available in NWChem's ESP module. You probably could write the code yourself, modifying ESP.

Sorry for the typo, of course not RESPA, but RESP.



I would like to do as described in the article, which introduced the RESP J Phys Chem 1993, v.97. p.10269
It describes three models with penalty functions:
(model 1) x2 = a (q₀ – q)²
(model 2) like (1) but q₀=0
(model 3) x2 = a ((q² +b²)^1/2 – b)

My first guess is that there is no control in NWChem for q₀. So I can use (model 2) and (model 3), but not (model 1). Is this true?

I see that the directive "constrain" does not provide the exact charges.
So my second assumption is that "constrain" not the charges fixing, but really is the (model 1). May be this true?

Okay, you may want to look for yourself at the code in esp_fit.F Looking at it, there is a definition for q0. Digging a little further and seeing where it is defined, turns out there is an not-advertized input variable called qh, which represents q0 in the RESP fitting. You can try setting:

  qh <real q0> 

in the esp block. I don't know why it is not included in the manual, may be it doesn't work properly. But, based on the equations you listed, it is the q0 you're looking for.

Bert

Thanks Bert.

I believe that

   qh <real q0>

is the total charge of the molecule or a charge for hydrogen.
Anyway this is not q0 separated per atoms.
So I'll write my own code. In this regards, one more question.
Do I understand correctly that the file nw.grid are the electrostatic potentials at the points on which the fit occur and only there, in the format

   x y z V

Looking at the code, it's a restraint on an atom (lines 235-265) but it's set the same for each atom, and it only is invoked for carbon and hydrogen. I would think this can easily be modified to suite your needs. If you develop a new functionality, you're more then welcome to contribute it back to the open-source development effort.

Yes, the grid file is x y z V.

Thanks,

Bert

Before you start coding it is key to understand the code, how it solves the equations you are interested in. You should not make guesses. What I recommend you do is to take the RESP paper and match up the coded equations to those in the paper. This will give you an understanding of what the code does and how it relates to the paper. Using the "q0" parameter should give you equation 1, at least when I match up the code. The b variable seems to match up with your b-variable in your equation 3.

I am not going to open files with binaries because of security issues. Have you done some tests that allow you to compare to the published data? This will allow you to test for correctness.

Bert

Unfortunately, I can not reproduce articles models 2 and 3 in NWChem too.
I have the impression that there is a problem with the units of “weight”.
Both the article and the documentation of NWChem says that it is the atomic units.
However, if you increase the hyperbolic constant in 50 times, and harmonic constant in 30 times, the results agree with the article.
Similarly, if the increase is 30 times the constant of my implementation (model 1), the results agree with the article.
So, I think that my (model 1) code is working properly, but it is necessary to deal with the input constants in all models.

The input file for methanol can be found here. Here are also modified sources. No binaries.
I allow the user defined target charges not only for harmonic functions, but for the hyperbolic too.

I didn't forget about you, was in China for a week and since I have been back I looked at this.

In general, I think your code should work fine. I included my updated and cleaned up version for you to test in File:Esp data.tgz.

It took me a long time to finally get a handle on the weight issue. If you look at the grid you define and the number of grid points you have a large number of grid points. However, in the paper you refer to, they only define a couple of radial shells, which end up at best to be about 2500 grid points for the potential to be fit on (the paper uses radial points at least 1.4 and at most 2.0 times the vdW radius). Why is this important, you should check the formulas in the paper. There is a dependency on the number of grid points. As a result, the more grid points you have the smaller the effect of the harmonic constant becomes.

I included an input deck with a different setup for defining the grid, which ends up closer to the on in the paper. We can never match the numbers in the paper, simply because we have a square cubic grid, while the paper uses spherical grids. However, with the grid points I generate in the input deck I do see the restraining effect as expected for the default constants.

Bert

Dear Bert,

Thank you very much for your help.
My task is solved now.


It seems to me that your update is not working. It falls on "esp" section, even if it is empty.

 Task  times  cpu:        0.5s     wall:        0.5s
0:0::: 0
(rank:0 hostname:abvgd pid:5936):ARMCI DASSERT fail. ../../ga-5-1/armci/src/common/armci.c:ARMCI_Error():208 cond:0


                                NWChem Input Module
                                -------------------


                     NWChem Electrostatic Potential Fit Module
                     -----------------------------------------
 geom_cart_get: geometry handle invalid            0
 geom_cart_get: open geometies:  1
  1 geom_cart_get: "geometry" -> "geometry"
 geom_cart_get: geometries in last accessed data base:  1
 geometry

 **********
 *   0: esp: geom_cart_get failed    0
 **********

If it has a matter to others, then this could be continued. If not, then I would stop there.