#14 new task

Vext implementation

Reported by: fgygi Owned by: fgygi
Priority: major Milestone:
Component: Qbox Version: 1.63.4
Keywords: vext Cc:


Implementation of the external potential interface.

Change History (8)

comment:1 Changed at 2016-08-22T21:14:12Z by fgygi

The change sets in patch file mahe-2016-08-17 have been checked in as r1874 (Basis.h) r1875 (qb.C) r1877 (RunCmd.C) and r1878 (all other files)

I have read the changes and I think most of the implementation is good. There are a number of points that will still require some modifications that I describe below. Please let me know what you think. We can address the following issues starting from r1878.

1) The use of vext as a pointer in Sample.h: I would like to keep this as a pointer. The simple interpretation is that if there is a non-null pointer, an external potential is used. Otherwise, if the pointer is null, there is no external potential. I prefer to not include an ExternalPotential member in Sample.h to avoid including additional dependencies on ExternalPotential.h in all classes that include Sample.h (as it is now, the only information needed is that ExternalPotential is a class, and the inclusion of ExternalPotential.h in Sample.h is not necessary.

2) Moving the vext file name in ExternalPotential: This is a good idea. Note that the current implementation of Qbox is not very consistent at this point regarding the location of the variables. Some values are in Control.h, and some are directly implemented as a member of some class. The change you proposed is more logical than it was before.

3) The use of the vext pointer in the Vext.h file: The proposed changes in Vext.h force the creation of an ExternalPotential object even if none is used. I prefer to avoid having resources allocated when they are not used, if possible (i.e. "you should not pay for what you don't need"). This means that we should have tests in other parts of the code (e.g. in EnergyFunctional.C) checking if the pointer s->vext is non-null before using the ExternalPotential object. Also, the print member of Vext would have to test the pointer before printing either s->vext->filename() or nothing.

4) ExternalPotential.h:
Using "using namespace std" statements in header files should be avoided (see e.g. discussions on issues of namespace protection on stackoverflow). The declarations of string and vector objects can be done using std::string and std::vector in header files. The "using namespace" statement is then used in the implementation file.

5) ExternalPotential.h: The use of reverse() to change the sign of the potential implies that the potential remains in a state where its sign depends on the history of use of the reverse() function. One possible way out of this is to use an amplitude variable when using the potential.

6) ResponseCmd.C: It is a good idea to implement the two different response calculations as two separate functions responseEfield and responseVext.
In the implementation of responseVext, it is not necessary to create new ChargeDensity objects (ChargeDensity construction is expensive, as it creates a number of other objects). The ChargeDensity member of BOSampleStepper can be accessed, and contains the information needed. We can save copies of the density in the variables rhor1 and rhor2 which would simply be vector<vector<double> > (which is cheaper than creating ChargeDensity objects).
At the end of the calculation of drho_r, the MPI_Barrier call should be limited to the communicator of the context (MPI_COMM_WORLD may affect other processes if the current context is smaller than COMM_WORLD). We will probably want to rewrite the communication of the data to the first-row nodes using a MPI_Gather call (and correspondingly a MPI_Scatter in the part where the data is read and distributed), but let's do that later.

7) RunCmd.C: Good catch, this had to be fixed.

Please add any comment in the comment section below.

comment:2 Changed at 2016-08-25T04:29:35Z by mahe

Dear Francois,

Thanks for your comments. I am looking at how to modify the code as you suggested.

Actually I have a long-standing question regarding the ChargeDensity class in Qbox. In my understanding, at any given time there should be only one ChargeDensity instance, why don't we create a charge density member in sample class? Logically a charge density should belong to a samply, just like wavefunction. And since Qbox is a DFT code where charge density is a critical component, I think it is reasonable to keep a charge density instance as long as Sample.wf exists.

When I am thinking about how to make ExternalPotential class more logical, I always found that the explicit dependence of update function on a charge density object to be not very natural.. If we could maintain a charge density instance during the whole process of simulation, then the ExternalPotential.update() could directly access the FourierTransform and Basis of charge density, which I think make sense..

I understand that you must have your considerations when implementing these and that it would be a lot of work if we change the way Qbox maintains charge density. I am just curious about the logic inside. Thanks!


comment:3 Changed at 2016-08-30T01:24:08Z by fgygi

I agree that it would be logical to have a ChargeDensity as a member of Sample. It would be interesting to consider the changes that would be needed to do that. As of now, it is only created by BOSampleStepper (or CPSampleStepper). The rationale for having it created "late" is that all members of Sample have to be updated whenever a command affects some of their characteristics (e.g. "set cell", "set ecut", "kpoint add" or "set nspin" would all trigger a "resize" operation on an existing ChargeDensity object. At the time it was written, it was avoided because of the cost (creating a new Basis, new FourierTransform at each kpoint, etc.). As an example, consider the scenario where an input file contains a long sequence of "kpoint add ..." commands. The ChargeDensity would have to be resized (including deleting and recreating FourierTransform objects after each "kpoint" command. (note that currently there is no resize() operation in ChargeDensity, since it is constructed with a fixed size that persists during its entire lifetime). It would be worth looking at the changes needed though, and ways to avoid costly resizes.

comment:4 Changed at 2016-08-30T01:37:59Z by fgygi

I have looked at the modifications of the vext branch. I checked in some modifications as r1879. Please have a look and let me know if you find there are problems and if you have another similar set of changes. This changeset seems to work with simple examples of external potentials I tried (linear potential on an H2O molecule in test/h2ogs).
At this point, the meaning of the vext variable is the following:

  • If vext is not defined, the pointer s->vext is null.
  • If vext is set, an ExternalPotential object is created and s->vext points to it.
  • If vext is set to NULL, the ExternalPotential is deleted.

The modifications do not resolve all issues. We still have to consider the following:
1) The distribution of the potential to the different tasks should be done with an MPI_Scatter operation rather than MPI_Bcast. Similarly, collecting the data of drho to task 0 could be simplified with MPI_Gather.
2) The semantics of "set vext filename" is not entirely clear. Should it update the external potential immediately? Should it only create an ExternalPotential object and initialize the filename_ member?
3) We have to consider what happens if the "response -vext v2.dat" command is used while there is already a vext defined using the "set vext v1.dat" command. In principle, this could be considered, meaning that the response to v2 on top of v1 would be computed. It makes the implementation rather complicated though, so we may want to keep that for later.

Please let me know about comments on r1879.

comment:5 Changed at 2016-08-30T01:47:14Z by fgygi

Two other issues we will have to consider:
4) Non-rectangular unit cells. This should not be a problem in principle, but there might be subtle points to consider.
5) We should define position of the origin of coordinates (in both Qbox and WEST) to have the same convention. In Qbox, the origin of coordinates of the grids does not coincide with an atom at position (0.0, 0.0, 0.0). That atom would be in the middle of the cell.

comment:6 Changed at 2016-09-05T03:27:14Z by mahe

Dear Francois,

Following are my notes for r1884 version.

  1. New functionality: 1.1. Multiple options added to response command. Now the way to call response command is:

response -vext vext_file [-RPA or -IPA] [-amplitude amplitude] [-parallel_write] nitscf [nite]
a) -RPA option will freeze the update of Vxc (only work for semilocal functional now); IPA will freeze the update of Vxc and Vhartree (not implemented yet since it will involve major change of EnergyFunctional class). In order to implment RPA, some changes are made to EnergyFunctional class and XCPotential class. Additionally, two boolean flags are added in Control class to indicate whether Vxc or Vhartree is freezed.
b) -amplitude option can set the amplitude of Vext. Right now it seems that this option is very convenient since it facilitates the test on whether the magnetitude for Vext is in linear-response regime. So I suggest keeping this option. c) -parallel_write: not implemented yet. By this option I intend to realize parallel I/O, like in SampleReader or SlaterDet::write. We may need to change the format of file into binary (e.g. base64) to implement this.

  1. Improvement of implementations: 2.1. Now response command will use the GS charge density and charge density under +Vext to generate an initial guess (by 2*rho0 - rho(+Vext)) for the calculation under -Vext. By doing this improvement it is observed that for silane the convergence of SCF iterations under -Vext becomes 3 times faster. In order to implement this there are some minor changes made to BOSamplerStepper. 2.2. Now in ExternalPotential::update process 0 will first read the file and MPI_Scatterv to other processes in column 0; then processes in column 0 will perform in-row MPI_Bcast to send Vext to all other processes. Likewise, now in ResponseCmd::responseVext process 0 will MPI_Gatherv drho_r from all processes on column 0. These functionalities are encapsulated in two functions in ExternalPotential.C and ResponseCmd.C respectively, to make the code more modular. As mentioned before, the IO performance might be further improved by resorting to binary format and MPI file operations.
  1. Other comments: 3.1. Regarding whether Vext should be updated when "set vext" command is called: right now our construction of Vext strongly depend on the ChargeDensity instance passed to update() function, so it seems that we can only update vext during "run" or "response" command. It would be interesting to consider if we can have a ChargeDensity member in Sample class, as you mentioned earlier. But it involves major changes. As a compromise, would it be possible if we keep a ChargeDensity member in Sample class (and we still have another ChargeDensity instance created in a BO sampler) and only update it manually by explicit command, or at the destruction of BO sampler? In this way, the awkwardness that Vext can only be updated inside a BO sampler where a charge density exists can be removed. 3.2. In ResponseCmd::responseVext function, a previous modification that changes s->wf.spincontext() to cd.context() seems to cause some tricky problem... Right now I rearlly don't understand why this will happen, but it appears that things work fine after I reverted this change. There might be some hidden bugs here that may be significant when we test against large calculations where the context is complicated. I will return to this problem later. 3.3. I understand that the option list of response command is too long compared to other commands in Qbox and it may be unpleasant. Later we may simplify it.


comment:7 Changed at 2016-09-07T22:25:01Z by mahe

A new revision is committed with some bugs fixed.


comment:8 Changed at 2016-09-11T00:34:56Z by fgygi

Updates r1895 and r1896 fix two errors affecting the response calculation.
r1896 reproduces the vext.dat.response of r1891 both in full and RPA mode.

Note: See TracTickets for help on using tickets.