Updated set_charges + set_initial_charges -stuff.

hotbit/aseinterface.py

@@ -375,8 +375,11 @@ def memory_estimate(self):
     def solve_ground_state(self,atoms):
         """ If atoms moved, solve electronic structure. """
         if not self.init:
+            assert type(atoms)!=type(None)
             self._initialize(atoms)
-        if self.calculation_required(atoms,'ground state'):
+        if type(atoms)==type(None):
+            pass
+        elif self.calculation_required(atoms,'ground state'):
             self.el.update_geometry(atoms)
             t0 = time()
             self.st.solve()
@@ -387,8 +390,8 @@ def solve_ground_state(self,atoms):
             self.flags['bonds'] = False
             if self.verbose:
                 print >> self.get_output(), "Solved in %0.2f seconds" % (t1-t0)
-            if self.get('SCC'):
-                atoms.set_charges(-self.st.get_dq())
+            #if self.get('SCC'):
+            #    atoms.set_charges(-self.st.get_dq())
         else:
             pass
 
@@ -790,14 +793,20 @@ def _init_mulliken(self):
             self.MA = MullikenAnalysis(self)
             self.flags['Mulliken']=True
 
-    def get_dq(self):
+    def get_dq(self,atoms=None):
         """ Return atoms' excess Mulliken populations.
         
         The total populations subtracted by
         the numbers of valence electrons.
         
         """
+        self.solve_ground_state(atoms)
         return self.st.get_dq()
+
+    def get_charges(self,atoms=None):
+        """ Return atoms' electric charges (Mulliken). """
+        return -self.get_dq(atoms)
+
 
     def get_atom_mulliken(self,I):
         """

hotbit/coulomb/baseclass.py

@@ -37,4 +37,5 @@ def get_forces(self, a=None):
         Return forces.
         """
         raise NotImplementedError()
-
+
+

hotbit/coulomb/direct_coulomb.py

@@ -58,7 +58,7 @@ def __init__(self, cutoff=None, timer=None):
 
     def update(self, a, q=None):
         if q is None:
-            q = a.get_charges()
+            q = a.get_initial_charges()
 
         r = a.get_positions()
         # FIXME!!! Check for change in cell, symmetries

hotbit/coulomb/multipole_expansion.py

@@ -107,7 +107,7 @@ def __init__(self, l_max=8, n=3, k=5, r0=None, timer=None):
 
     def update(self, a, q=None):
         if q is None:
-            q = a.get_charges()
+            q = a.get_initial_charges()
 
         r = a.get_positions()
         # FIXME!!! Check for change in cell, symmetries

hotbit/test/madelung_constants.py

@@ -75,18 +75,18 @@
 
         #a.set_charges([ (Q if sym == syms[0] else -Q) for sym in syms ])
         # to work with older versions
-        a.set_charges([ (-Q,Q)[sym==syms[0]] for sym in syms ])
+        a.set_initial_charges([ (-Q,Q)[sym==syms[0]] for sym in syms ])
 
         a.translate([0.25*a0,0.25*a0,0.25*a0])
         if debug:
             write("%s.cfg" % name, a)
 
         ha = HotbitAtoms(a, container='Bravais')
 
-        sol.update(ha, ha.get_charges())
+        sol.update(ha, ha.get_initial_charges())
 
         phi  = sol.get_potential()
-        e    = np.sum(a.get_charges()*phi)/2
+        e    = np.sum(a.get_initial_charges()*phi)/2
         M    = -2*e*a0*nnd/(len(a))
         err  = abs(M-target_M)
 

hotbit/test/misc.py

@@ -1,5 +1,5 @@
 import os
-from ase.data.molecules import molecule as ase_molecule
+from ase.structure import molecule as ase_molecule
 from ase import *
 from hotbit import fixpar
 

hotbit/test/multipole_operations.py

@@ -71,7 +71,7 @@
     # Compute moments
     M = [ ]
     for i in range(8):
-        M += [ get_moments(a[i].get_positions(), a[i].get_charges(), L_MAX, r0) ]
+        M += [ get_moments(a[i].get_positions(), a[i].get_initial_charges(), L_MAX, r0) ]
 
     # Construct a composite atoms object
     # and compute the corresponding multipole
@@ -93,7 +93,7 @@
         multipole_to_multipole(dr, L_MAX, M0_l, M_L, Mc0_l, Mc_L)
 
     # Compute the multipole moment directly
-    Md0_l, Md_L = get_moments(b.get_positions(), b.get_charges(), L_MAX, r0c)
+    Md0_l, Md_L = get_moments(b.get_positions(), b.get_initial_charges(), L_MAX, r0c)
 
     err_mom1 = np.max(np.abs(Mc0_l-Md0_l))
     err_mom2 = np.max(np.abs(Mc_L-Md_L))
@@ -158,7 +158,7 @@
     assert err_phi4 < TOL_PHI2
 
     # Compute the multipole moment directly
-    Md0_l, Md_L = get_moments(b.get_positions(), b.get_charges(), L_MAX, r0c)
+    Md0_l, Md_L = get_moments(b.get_positions(), b.get_initial_charges(), L_MAX, r0c)
 
     # Now rotate the atoms, recompute the multipole expansion and compare the result
     # to the rotated expansion
@@ -180,7 +180,7 @@
     for i in b:
         i.position = r0c + np.dot(R, i.position - r0c)
 
-    Me0_l, Me_L = get_moments(b.get_positions(), b.get_charges(), L_MAX, r0c)
+    Me0_l, Me_L = get_moments(b.get_positions(), b.get_initial_charges(), L_MAX, r0c)
     Mf0_l, Mf_L = transform_multipole(R, L_MAX, Md0_l, Md_L)
 
     err_mom3 = np.max(np.abs(Me0_l-Mf0_l))

hotbit/test/periodicity.py

@@ -52,7 +52,7 @@ def electrostatics_test(b, r=3, r0=None):
 
     # Check multipole moments
     mp  = MultipoleExpansion(L_MAX, r, k3, r0)
-    mp.update(b, b.get_charges())
+    mp.update(b, b.get_initial_charges())
 
     # Store moments and field for later
     moments          = mp.get_moments()
@@ -86,7 +86,7 @@ def electrostatics_test(b, r=3, r0=None):
 
     # Next neighbor shell by direct summation
     mp  = MultipoleExpansion(L_MAX, r*r*r, k1, r0)
-    mp.update(b, b.get_charges())
+    mp.update(b, b.get_initial_charges())
 
     phi_dir, E_dir  = mp.get_potential_and_field(b)
 
@@ -98,7 +98,7 @@ def electrostatics_test(b, r=3, r0=None):
     rep  = [ ((0,0),(-(r-1)/2, (r-1)/2))[i] for i in b.get_pbc() ]
     c    = b.extended_copy(tuple(rep))
 
-    M0_l, M_L  = get_moments(c.get_positions(), c.get_charges(), L_MAX, r0)
+    M0_l, M_L  = get_moments(c.get_positions(), c.get_initial_charges(), L_MAX, r0)
 
     #print M_L
     #print M_L_mp

hotbit/test/polyethene_twisted.py

@@ -59,7 +59,7 @@
     atoms.rattle(0.1)
 
     # Check electrostatics only
-#    atoms.set_charges([1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0])
+#    atoms.set_initial_charges([1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0])
 #    atoms.set_calculator(calc.st.es)
 
     # Check forces from finite differences

hotbit/test/standard_set.py

@@ -1,7 +1,7 @@
 from hotbit import Hotbit 
 #from hotbit import Calculator0
 from ase import *
-from ase.data.molecules import molecule
+from ase.structure import molecule
 from hotbit.test.misc import default_param
 import sys