esp restrain


Gets Around
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?

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

Bert


Quote:P99 Oct 8th 11:53 am
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?

Gets Around
I would not want to fix the charges exactly, but to use RESPA.

Quote:Bert Oct 8th 9:31 am
The options in constrains (just above restraints in the documentation) should be able to allow you to do what you want I believe.

Bert

Forum Vet
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.

Bert

[QUOTE=P99 Oct 8th 5:02 pm]I would not want to fix the charges exactly, but to use RESPA.

Quote:Bert Oct 8th 9:31 am
The options in constrains (just above restraints in the documentation) should be able to allow you to do what you want I believe.

Bert

Gets Around
Quote:Bert Oct 9th 1:16 pm
I have no clue what RESPA stands for and can't find anything online about it either.

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

Quote:Bert Oct 9th 1:16 pm
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.


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 (q0 – q)2
(model 2) like (1) but q0=0
(model 3) x2 = a ((q2 +b2)1/2 – b)

My first guess is that there is no control in NWChem for q0. 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?

Forum Vet
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


Quote:P99 Oct 11th 9:37 am
Quote:Bert Oct 9th 1:16 pm
I have no clue what RESPA stands for and can't find anything online about it either.

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

Quote:Bert Oct 9th 1:16 pm
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.


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 (q0 – q)2
(model 2) like (1) but q0=0
(model 3) x2 = a ((q2 +b2)1/2 – b)

My first guess is that there is no control in NWChem for q0. 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?

Gets Around
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

Quote:Bert Oct 11th 12:48 pm
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

Forum Vet
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



Quote:P99 Oct 12th 8:49 am
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

Quote:Bert Oct 11th 12:48 pm
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

Gets Around
Hi

It seems to me that I have implemented (method 1). File with the modified source, nwchem.exe and example linked.

By analogy with the "radius" I introduced in "ESP" module a command "q0":

      q0 <integer iatnum> <real q0>
where iatnum is the atomic number for which a q0 will be used.

However, there are some doubts.

1. I just followed the analogy without real understanding what is written in the code.

First, I do not understand what an undocumented variable "q0" do.
It does something strange... plus activates RESP3, whatever that meant.
So I did not used it, but introduced a new variable "q00" (array).

Secondly, I do not understand what is the b(j) variable. I simply ignored it.

(My changes are marked in the text by commentary "NIKI")

2. I do not know how global variables work in Fortran. So I made a file niki.fh. Although it looks like working, however, it is an obvious hack.


Could you be so kind as to look at my source code?

(Sorry I am working under Windows, so executable is a Windows exe file.)

Quote:Bert Oct 12th 12:09 pm
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


Forum Vet
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

Gets Around
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.



Quote:Bert Oct 17th 7:44 pm
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

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


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

 

and
Quote:P99 Oct 18th 12:24 pm
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.



Quote:Bert Oct 17th 7:44 pm
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

Gets Around
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.



Quote:Bert Oct 26th 3:02 pm
I didn't forget about you, was in China for a week and since I have been back I looked at this.


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

Forum Vet
This is related to the different versions of the code we're using I believe. So, this issue is closed.

Bert


Quote:P99 Oct 29th 2:05 pm
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.



Quote:Bert Oct 26th 3:02 pm
I didn't forget about you, was in China for a week and since I have been back I looked at this.


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

Gets Around
Hi,

I compiled Bert update with the current version of NWChem (6.3). With it, it works and it seems to me that it works well. The only disadvantage is that the 'q0' option became mandatory.

Could you be so kind to enter this code in the standard NWChem?
It seems to me that the presence of model 1 of J Phys Chem 1993, v.97. p.10269 a useful addition.


Forum >> NWChem's corner >> General Topics