[ase-users] VASP module in ASE

Jens Jørgen Mortensen jensj at fysik.dtu.dk
Mon Sep 1 12:53:21 CEST 2008


On Fri, 2008-08-29 at 08:01 +0200, Jonas Björk wrote:
> Hi,
> 
> I have come up with a new version of the VASP module by bringing
> together Jussi Enkovaara's and John Kitchin's code and performing some
> modifications and adding some extra features myself. 

Thanks.  I have a few comments - see below ...

> I have made some changes as I discussed with Jens Jørgen:
> For example if a non-zero magnetic moment is specified for any atom a
> spin-polarized calculation will be performed. One thing I have marked
> out in the calculator is the "get_magnetic_moments" because the
> calculated magnetic moments for each atom in vasp doesn't add up to the
> total magnetic moment because it is defined as the magnetic moment
> inside a wigner-seitz radius in some way ... Furthermore to get the
> magnetic moments for each atoms one needs so specify some extra
> parameters in the input file, which gives extra output files which are
> totally unnecessary for the purpose. I could turn on this feature again
> but then it would only return the magnetic moment for the individual
> atoms if the user have specified the necessary parameter(s) himself. 
> The total magnetic moment is calculated fine though. 

I think it's good to have the get_magnetic_moments() method.  We can
just write in its docstring which VASP parameters to set in order to
make it work.

> I have kept the different input defined in 'incar_parameters' from
> Jussi's code, which defines the parameters which should be written to
> the vasp input (INCAR) file (except the parameter for initial magnetic
> moment and spin-polarized calculation which is decided from the atoms
> object). Other parameters are saved in 'input_parameters'. 
> 
> A parameter can be set either when the calculator is called, e.g.
> 
> calc=Vasp(nbands=24)
> 
> or as
> 
> calc.set(nbands=24)
> 
> Should there be a special command for specifying the kpts or is
> 
> calc=Vasp(kpts=...), calc.set(kpts=...) alright?

That's fine.

> The same for the
> planewave cutoff. In Vasp it is called 'encut'. Should there be a
> function like calc.set_planewave_cutoff() or something?

No.  Let's just use 'encut'.

> Regarding the k-points they can either be specified by the number of
> k-points, or explicitly k-point per k-point. The interface decide
> automatically which action to take when writing the Vasp KPOINTS file. 

Perfect!

> The POTCAR file can now be built either using LDA, PW91 or PBE xc, I
> will also add the support for the ultrasoft PPs.
> 
> Please come with commands on what should be changed, added or taken away
> from the code. I have added quite a few functions (e.g. get_eigenvalues)
> which not are implemented yet, most of them will come I hope. 

That's fine - let's get the basic stuff working first.

I made some changes to the code (vasp3.py).  Below is a diff with some
comments marked with ???.  If you send me a corrected version - then I
will check it in.

Jens Jørgen

[jensj at casimir ~]$ diff -u vasp2.py vasp3.py
--- vasp2.py    2008-09-01 12:18:55.778666000 +0200
+++ vasp3.py    2008-09-01 12:36:53.688641000 +0200
@@ -16,7 +16,7 @@
 
 
 class Vasp:
-    def __init__(self, restart=False, **kwargs):
+    def __init__(self, restart=None, **kwargs):

??? use None for default values.

         # Parameters that can be set in INCAR. The values which are
None
         # are not written and default parameters of VASP are used for
them.
@@ -98,15 +98,17 @@
 
         if restart:
             self.read_incar()
-            self.atoms=ase.io.read('CONTCAR', format='vasp')
+            self.atoms = ase.io.read('CONTCAR', format='vasp')

??? Add space arround '=' for readability and to conform with our ASE
standard:

https://wiki.fysik.dtu.dk/ase/development/code.html#python-coding-conventions

http://www.python.org/dev/peps/pep-0008/

 #            self.read_kpoints()
             return
 
         if self.input_parameters['xc'] not in ['91','lda','pbe']:
-            raise Exception, "%s not supported for xc! use one of
['91',lda','pbe']" % kwargs['xc']
-        self.positions=None
-        self.nbands=self.incar_parameters['nbands']
-        self.atoms=None
+            raise Exception(
+                "%s not supported for xc! use one of ['91',lda','pbe']"
%
+                kwargs['xc'])

??? use the syntax Exception('text')

+        self.positions = None
+        self.nbands = self.incar_parameters['nbands']
+        self.atoms = None
         self.set(**kwargs)
 
 
@@ -117,9 +119,7 @@
             elif self.incar_parameters.has_key(key):
                 self.incar_parameters[key]=kwargs[key]
             else:
-                print 'Parameter '+key+' not defined'
-                return
-
+                raise TypeError('Parameter not defined: ' + key)

??? A wrong key should raise an exception.
 
     def update(self, atoms):
         if (self.positions is None or
@@ -132,7 +132,6 @@
             self.initialize(atoms)
             self.calculate(atoms)
 
-

??? Use just a single blank line to separate methods.

     def initialize(self, atoms):
         """Initialize a VASP calculation
 
@@ -146,22 +145,25 @@
         self.spinpol=atoms.get_initial_magnetic_moments().any()
         # Determine the number of atoms of each atomic species
         num = atoms.get_chemical_symbols()
-        atom_num=[[num[0],1]]
-        for m in range(1,len(num)):
-            if num[m]==atom_num[-1][0]:
-                atom_num[-1][1]+=1
-            else:
-                atom_num.append([num[m], 1])
+
+        # How can this ever work?  In ASE, we don't assume that the
+        # atoms are ordered in any way?
+        atom_num=[[num[0],1]]                 ???
+        for m in range(1,len(num)):           ???
+            if num[m]==atom_num[-1][0]:       ???
+                atom_num[-1][1]+=1            ???
+            else:                             ???
+                atom_num.append([num[m], 1])  ???

??? Here you are assuming that the atoms are ordered after atomic number
- right?  We should fix that!

         self.atomlist=[]
         self.symbols=[]
         for n in range(len(atom_num)):
             self.atomlist.append(atom_num[n][1])
             self.symbols.append(atom_num[n][0])
         xc='/'
-        if p['xc']=='91':
-            xc='_gga/'
-        elif p['xc']=='_pbe/':
-            xc=pbe
+        if p['xc'] == 'PW91':
+            xc = '_gga/'
+        elif p['xc'] == 'PBE':
+            xc = '_pbe/'

??? How about using the the names LDA, PW91, and PBE?

         if 'VASP_PP_PATH' in os.environ:
             pppaths = os.environ['VASP_PP_PATH'].split(':')
         else:
@@ -191,10 +193,9 @@
         self.converged = False
         self.setups_changed = False
 
-
     def calculate(self, atoms):
-        """
-        Method that generate necessary files in the working directory. 
+        """Generate necessary files in the working directory.
+        
         If the directory does not exist it will be created.
         """

??? A docstring must start with stand-alone headline - possibly followed
by a longer description.

         positions = atoms.get_positions()
@@ -218,25 +219,23 @@
         if exitcode != 0:
             raise RuntimeError('Vasp exited with exit code: %d.  ' %
exitcode)
         
-        self.positions=positions.copy(atoms)
-        [self.energy_free, self.energy_zero]=self.read_energy()
-        self.forces=self.read_forces(atoms)
+        self.positions = positions.copy(atoms)
+        self.energy_free, self.energy_zero = self.read_energy()
+        self.forces = self.read_forces(atoms)

??? More cosmetics

         self.dipole=self.read_dipole()
         self.fermi=self.read_fermi()
         self.atoms=atoms.copy()
         if not self.nbands:
             self.nbands=self.read_nbands()
         if self.spinpol:
-            self._magnetic_moment=self.read_magnetic_moment()
+            self.magnetic_moment = self.read_magnetic_moment()

??? Don't use _names_starting_with_underscores unless there is a very
good reason for it.

         self.set(nbands=self.nbands)
         self.old_incar_parameters=self.incar_parameters.copy()
         self.old_input_parameters=self.input_parameters.copy()
 
-
     def get_atoms(self):
         return self.atoms
 
-
     def get_potential_energy(self, atoms, force_consistent=False):
         self.update(atoms)
         if force_consistent:
@@ -244,16 +243,13 @@
         else:
             return self.energy_zero
 
-
     def get_forces(self, atoms):
         self.update(atoms)
         return self.forces
 
 
     def get_stress(self, atoms):
-        #Not implemented yet
-        return
-
+        raise NotImplementedError

??? Use the NotImplementedError exception for such cases

     def calculation_required(self,atoms, quantities):
         # Not implemented yet
@@ -401,8 +397,7 @@
 
 
 
-    """ Read information from OUTCAR files"""
-
+    # Methods for reading information from OUTCAR files:
     def read_energy(self):
         file=open('OUTCAR','r')
         lines=file.readlines()


-------------- next part --------------
A non-text attachment was scrubbed...
Name: vasp3.py
Type: text/x-python
Size: 17647 bytes
Desc: not available
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20080901/087e152a/attachment.py>


More information about the ase-users mailing list