implemented multiple ncsf VASP cycles

doc/tutorials/images_scripts/nio_csc.py

@@ -105,7 +105,7 @@ def dmft_cycle():
         ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
         ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.cthyb_hash
         ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
-        ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.dft_tools_hash
+        ar['DMFT_input']['code_versions']["dft_tools_git"] = dft_tools_version.dft_tools_hash
         if 'iteration_count' in ar['DMFT_results']: 
             iteration_offset = ar['DMFT_results']['iteration_count']+1
             S.Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]

doc/tutorials/nio_csc.rst

@@ -53,7 +53,7 @@ DMFT
 dmft script
 ------------------------------
 
-Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on :ref:`srvo3`, we will not go into detail here but simply provide the script :ref:`nio.py`. Following Kunes et al. in `PRB 75 165115 (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.165115>`_ we use :math:`U=8` and :math:`J=1`. We se;ect :math:`\beta=5` instead of :math:`\beta=10` to ease the problem slightly. For simplicity we fix the double-counting potential to :math:`\mu_{DC}=59` eV by::
+Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on :ref:`SrVO3 <srvo3>`, we will not go into detail here but simply provide the script :ref:`nio.py`. Following Kunes et al. in `PRB 75 165115 (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.165115>`_ we use :math:`U=8` and :math:`J=1`. We se;ect :math:`\beta=5` instead of :math:`\beta=10` to ease the problem slightly. For simplicity we fix the double-counting potential to :math:`\mu_{DC}=59` eV by::
 
   DC_value = 59.0
   SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_value=DC_value)
@@ -62,11 +62,11 @@ For sensible results run this script in parallel on at least 20 cores. As a quic
 
 
 Local lattice Green's function for all projected orbitals
-----------------------
+---------------------------------------------------------
 We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script :ref:`NiO_local_lattice_GF.py`. The result is saved in the h5 file as `G_latt_orb_it<n_it>`, where `<n_it>` is the number of the last DMFT iteration.
 
 Spectral function on real axis: MaxEnt
-----------------------
+--------------------------------------
 To compare to results from literature we make use of the `maxent triqs application <https://triqs.github.io/maxent/master/>`_ and calculate the spectral function on real axis. Use this script to perform a crude but quick calculation: :ref:`maxent.py` using a linear real axis and a line-fit analyzer to determine the optimal :math:`\alpha`. The output is saved in the h5 file in `DMFT_results/Iterations/G_latt_orb_w_o<n_o>_it<n_it>`, where `<n_o>` is the number of the orbital and `n_it` is again the number of the last iteration. The real axis information is stored in `DMFT_results/Iterations/w_it<n_it>`.
 
 
@@ -74,26 +74,28 @@ To compare to results from literature we make use of the `maxent triqs applicati
     :width: 400
     :align: center
 
-.. Charge self-consistent DMFT
+Charge self-consistent DMFT
-.. ==================================================
+==================================================
 
 
-.. In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP `INCAR` such that :program:`VASP` reads the updated charge density after each step. We add the lines
+In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP `INCAR` such that :program:`VASP` reads the updated charge density after each step. We add the lines::
 
-..   ICHARG = 5
+  ICHARG = 5
-..   NELM = 1000
+  NELM = 1000
-..   NELMIN = 1000
+  NELMIN = 1000
 
-.. which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file `GAMMA`. It is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number.
+which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file `GAMMA`. It is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number.
 
-.. We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `CHGCAR` and `nio.h5` together with the other :program:`VASP` input files and :file:`plo.cfg` in a new directory. We use a script called :program:`vasp_dmft` to invoke :program:`VASP` in the background and start the DMFT calculation together with :program:`plovasp` and the converter. This script assumes that the dmft sript contains a function `dmft_cycle()` and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforardly and call it :ref:`nio_csc.py`. The most important new lines are
+We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `CHGCAR` and `nio.h5` together with the other :program:`VASP` input files and :file:`plo.cfg` in a new directory. We use a script called :program:`vasp_dmft` to invoke :program:`VASP` in the background and start the DMFT calculation together with :program:`plovasp` and the converter. This script assumes that the dmft sript contains a function `dmft_cycle()` and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforardly (for a working example see :ref:`nio_csc.py`). The most important new lines are::
 
-..  SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
+  SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
-..  SK.calc_density_correction(dm_type='vasp')
+  SK.calc_density_correction(dm_type='vasp')
-..  correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
+  correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
 
-.. where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file :file:`GAMMA`. The correlation energy is calculated via Migdal-Galitzki formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.
+where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file :file:`GAMMA`. The correlation energy is calculated via Migdal-Galitzki formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.
 
-.. We can start the whole machinery by excectuing
+To help convergence, we keep the density (i.e., the GAMMA file) fixed for a few DFT iterations. We do so since VASP uses an iterative diagonalization.
 
-..  vasp_dmft -n <n_procs> -i <n_iters> -p <vasp_exec>  nio_csc.py
+We can start the whole machinery by excecuting::
+
+  vasp_dmft -n <n_procs> -i <n_iters> -j <n_iters_dft> -p <vasp_exec>  nio_csc.py

python/converters/plovasp/elstruct.py

@@ -140,6 +140,7 @@ class ElectronicStructure:
         ions = sorted(list(set([param['isite'] for param in self.proj_params])))
         nions = len(ions)
         norb = nproj / nions
+        print nproj, ns, nk, nb, ions, nions, norb
 
 # Spin factor
         sp_fac = 2.0 if ns == 1 and not self.nc_flag else 1.0
@@ -162,12 +163,15 @@ class ElectronicStructure:
         for ispin in xrange(ns):
             print
             print "  Spin:", ispin + 1
+            print self.proj_params
             for io, ion in enumerate(ions):
                 print "  Site:", ion
                 iorb_inds = [(ip, param['m']) for ip, param in enumerate(self.proj_params) if param['isite'] == ion]
+                print iorb_inds
                 norb = len(iorb_inds)
-                dm = np.zeros((norb, norb))
+                norb2 = self.proj_params[0]['l']*2+1
-                ov = np.zeros((norb, norb))
+                dm = np.zeros((norb2, norb2))
+                ov = np.zeros((norb2, norb2))
                 for ind, iorb in iorb_inds:
                     for ind2, iorb2 in iorb_inds:
                         dm[iorb, iorb2] = den_mat[ispin, ind, ind2]

python/converters/plovasp/sc_dmft.py

@@ -32,6 +32,7 @@ import signal
 import sys
 import pytriqs.utility.mpi as mpi
 import converter
+from shutil import copyfile
 
 xch = sys.excepthook
 def excepthook(typ, value, traceback):
@@ -40,7 +41,7 @@ def excepthook(typ, value, traceback):
         mpi.MPI.COMM_WORLD.Abort(1)
 sys.excepthook = excepthook
 
-debug = True
+debug = False
 #
 # Helper functions
 #
@@ -99,7 +100,7 @@ class bcolors:
 #
 # Main self-consistent cycle
 #
-def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
+def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
     """
     """
     mpi.report("  Waiting for VASP lock to appear...")
@@ -120,6 +121,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
                 mpi.report("  VASP stopped")
                 vasp_running = False
                 break
+            
 # Tell VASP to stop if the maximum number of iterations is reached
         iter += 1
         if iter == n_iter:
@@ -155,8 +157,29 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
             print "="*80
             print
         
-        if mpi.is_master_node() and vasp_running:
+        # check if we should do additional VASP calculations
+        # in the standard VASP version, VASP writes out GAMMA itself
+        # so that if we want to keep GAMMA fixed we have to copy it to
+        # GAMMA_recent and copy it back after VASP has completed an iteration
+        # if we are using a hacked Version of VASP where the write out is
+        # disabled we can skip this step.
+        # the hack consists of removing the call of LPRJ_LDApU in VASP src file
+        # electron.F around line 644
+        iter_dft = 0
+        if vasp_version == 'standard':
+            copyfile(src='GAMMA',dst='GAMMA_recent')
+        while iter_dft < n_iter_dft:
+            if mpi.is_master_node():
                 open('./vasp.lock', 'a').close()
+            while is_vasp_lock_present():
+                time.sleep(1)
+                if not is_vasp_running(vasp_pid):
+                    mpi.report("  VASP stopped")
+                    vasp_running = False
+                    break
+            iter_dft += 1
+            if vasp_version == 'standard':
+                copyfile(src='GAMMA_recent',dst='GAMMA')
                     
     if mpi.is_master_node():
         total_energy = dft_energy + corr_energy - dft_dc
@@ -170,7 +193,9 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
     mpi.report("***Done")
 
 def main():
+    
     import importlib
+    
     try:
         vasp_pid = int(sys.argv[1])
     except (ValueError, KeyError):
@@ -186,18 +211,33 @@ def main():
         raise
         
     try:
-        dmft_script = re.sub("\.py$", "", sys.argv[3])
+        n_iter_dft = int(sys.argv[3])
+    except (ValueError, KeyError):
+        if mpi.is_master_node():
+            print "ERROR: Number of VASP iterations with fixed charge density must be provided as the third argument"
+        raise
+
+    try:
+        dmft_script = re.sub("\.py$", "", sys.argv[4])
     except:
         if mpi.is_master_node():
-            print "ERROR: User-defined DMFT script must be provided as the third argument"
+            print "ERROR: User-defined DMFT script must be provided as the fourth argument"
         raise
 
 # Optional parameter: config-file name
     try:
-        cfg_file = sys.argv[4]
+        cfg_file = sys.argv[5]
     except KeyError:
         cfg_file = 'plo.cfg'
         
+    try:
+        vasp_version = sys.argv[6]
+    except KeyError:
+        vasp_version = 'standard'
+    
+    if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
+        raise Exception('vasp_version has to be standard or no_gamma_write')
+
 #    if len(sys.argv) > 1:
 #        vasp_path = sys.argv[1]
 #    else:
@@ -211,7 +251,7 @@ def main():
 
     dmft_mod = importlib.import_module(dmft_script)
 
-    run_all(vasp_pid, dmft_mod.dmft_cycle, cfg_file, n_iter)
+    run_all(vasp_pid, dmft_mod.dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version)
 
 if __name__ == '__main__':
     main()

shells/vasp_dmft.bash.in

@@ -5,9 +5,9 @@ MPIRUN_CMD=mpirun
 show_help()
 {
 echo "
-Usage: vasp_dmft [-n <number of cores>] -i <number of iterations> [-p <path to VASP directory>] [<dmft_script.py>]
+Usage: vasp_dmft [-n <number of cores>] -i <number of iterations>  -j <number of VASP iterations with fixed charge density> [-v <VASP version>] [-p <path to VASP directory>] [<dmft_script.py>]
 
-       If the number of cores is note specified it is set to 1 by default.
+       If the number of cores is not specified it is set to 1 by default.
 
        <dmft_script.py> must provide an importable function 'dmft_cycle()' which is invoked
        once per DFT+DMFT iteration. If the script name is omitted the default name 
@@ -18,7 +18,7 @@ Usage: vasp_dmft [-n <number of cores>] -i <number of iterations> [-p <path to V
 "
 }
 
-while getopts ":n:i:p:" opt; do
+while getopts ":n:i:j:v:p:h" opt; do
   case $opt in
      n)
 #        echo "Option: Ncpu, argument: $OPTARG"
@@ -34,11 +34,27 @@ while getopts ":n:i:p:" opt; do
 #           echo "Number of iterations: $NITER"
         fi
         ;;
+     j)
+#        echo "Option: Ndftiter"
+        if [ -n "$OPTARG" ]; then
+           NDFTITER=$OPTARG
+#           echo "Number of iterations with fixed density: $NDFTITER"
+        fi
+        ;;
      p)
         if [ -n "$OPTARG" ]; then
            VASP_DIR=$OPTARG
         fi
         ;;
+     v)
+        if [ -n "$OPTARG" ]; then
+           VASP_VERSION=$OPTARG
+        fi
+        ;;
+      h)
+        show_help
+        exit 1
+        ;;
       :)
         echo "  Error: Option -$OPTARG requires an argument" >&2
         show_help
@@ -55,6 +71,7 @@ if [ -z "$NITER" ]; then
   exit 1
 fi
 
+
 if [ -z "$VASP_DIR" ]; then
   echo "  Error: A path to VASP directory must be given either with option -p or by setting variable VASP_DIR" >&2
   show_help
@@ -66,6 +83,16 @@ if [ -z "$NPROC" ]; then
   NPROC=1
 fi
 
+if [ -z "$NDFTITER" ]; then
+  echo "  Number of VASP iterations without updating density not specified, setting to 1"
+  NDFTITER=1
+fi
+
+if [ -z "$VASP_VERSION" ]; then
+  echo "  VASP version not specified, setting to standard"
+  VASP_VERSION="standard"
+fi
+
 shift $((OPTIND-1))
 
 if [ -z "$1" ]; then
@@ -76,11 +103,14 @@ fi
 
 echo "  Number of cores: $NPROC"
 echo "  Number of iterations: $NITER"
+echo "  Number of iterations with fixed density: $NDFTITER"
+echo "  VASP version: $VASP_VERSION"
 echo "  Script name: $DMFT_SCRIPT"
 
 rm -f vasp.lock
-stdbuf -o 0 $MPIRUN_CMD -np $NPROC "$VASP_DIR" &
+# stdbuf -o 0 $MPIRUN_CMD -np $NPROC "$VASP_DIR" &
+vasp_intel.sh &
 
 
-$MPIRUN_CMD -np $NPROC python -m triqs_dft_tools.converters.plovasp.sc_dmft $(jobs -p) $NITER $DMFT_SCRIPT 'plo.cfg' || kill %1
+$MPIRUN_CMD -np $NPROC python -m triqs_dft_tools.converters.plovasp.sc_dmft $(jobs -p) $NITER $NDFTITER $DMFT_SCRIPT 'plo.cfg' $VASP_VERSION || kill %1