[ase-users] NPTBerendsen()
Vitaly Vitalievich
vvchaban at gmail.com
Mon Oct 28 10:38:06 CET 2024
I hereby share certain modifications I made to the MOPAC() class should
anyone volunteer to implement the NPT functionality.
def read_results(self):
"""Read the results, such as energy, forces, eigenvalues, etc.
"""
FileIOCalculator.read(self, self.label)
if not os.path.isfile(self.label + '.out'):
raise ReadError
with open(self.label + '.out') as fd:
lines = fd.readlines()
for i, line in enumerate(lines):
if line.find('TOTAL ENERGY') != -1:
self.results['energy'] = float(line.split()[3])
elif line.find('FINAL HEAT OF FORMATION') != -1:
self.final_hof = float(line.split()[5]) * kcal / mol
self.results['energy'] = self.final_hof
elif line.find('NO. OF FILLED LEVELS') != -1:
self.nspins = 1
self.no_occ_levels = int(line.split()[-1])
elif line.find('NO. OF ALPHA ELECTRON') != -1:
self.nspins = 2
self.no_alpha_electrons = int(line.split()[-1])
self.no_beta_electrons = int(lines[i + 1].split()[-1])
self.results['magmom'] = abs(self.no_alpha_electrons -
self.no_beta_electrons)
elif line.find('FINAL POINT AND DERIVATIVES') != -1:
forces = [-float(line.split()[6])
for line in lines[i + 3:i + 3 + 3 *
len(self.atoms)]]
self.results['forces'] = np.array(
forces).reshape((-1, 3)) * kcal / mol
forces = [-float(line.split()[6])
for line in lines[i + 3:i + 3 + 3 *
len(self.atoms)]]
self.results['forces'] = np.array(
forces).reshape((-1, 3)) * kcal / mol
* elif line.find('VOLUME OF UNIT CELL') != -1:
self.volume = float(line.split()[5]) * units.Ang^3
self.results['volume'] = self.volume *
* elif line.find('Stress tensor in GPa using Voigt notation') !=
-1: if self.atoms.cell.rank == 3:
self.stress += self.stress.T.copy() self.stress *= -0.5
/ self.atoms.get_volume() self.results['stress'] =
self.stress.flat[[0, 4, 8, 5, 2, 1]]*
elif line.find('EIGENVALUES') != -1:
if line.find('ALPHA') != -1:
j = i + 1
eigs_alpha = []
while not lines[j].isspace():
eigs_alpha += [float(eps) for eps in
lines[j].split()]
j += 1
elif line.find('BETA') != -1:
j = i + 1
eigs_beta = []
while not lines[j].isspace():
eigs_beta += [float(eps) for eps in
lines[j].split()]
j += 1
eigs = np.array([eigs_alpha, eigs_beta]).reshape(2, 1,
-1)
self.eigenvalues = eigs
else:
eigs = []
j = i + 1
while not lines[j].isspace():
eigs += [float(e) for e in lines[j].split()]
j += 1
self.eigenvalues = np.array(eigs).reshape(1, 1, -1)
elif line.find('DIPOLE ') != -1:
self.results['dipole'] = np.array(
lines[i + 3].split()[1:1 + 3], float) * Debye
def get_eigenvalues(self, kpt=0, spin=0):
return self.eigenvalues[spin, kpt]
def get_homo_lumo_levels(self):
eigs = self.eigenvalues
if self.nspins == 1:
nocc = self.no_occ_levels
return np.array([eigs[0, 0, nocc - 1], eigs[0, 0, nocc]])
else:
na = self.no_alpha_electrons
nb = self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
eah, eal = eigs[0, 0, na - 1: na + 1]
ebh, ebl = eigs[1, 0, nb - 1: nb + 1]
return np.array([max(eah, ebh), min(eal, ebl)])
def get_somo_levels(self):
assert self.nspins == 2
na, nb = self.no_alpha_electrons, self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
return np.array([self.eigenvalues[0, 0, na - 1],
self.eigenvalues[1, 0, nb - 1]])
def get_final_heat_of_formation(self):
return self.final_hof
* def get_volume (self): return self.volume*
On Mon, Oct 28, 2024 at 11:25 AM Vitaly Vitalievich <vvchaban at gmail.com>
wrote:
> Adam, thank you for your insights.
>
> Here is how the stress looks in MOPAC22. If you happen to know the
> internals of ASE deeply, unlike me, it must be straightforward to code.
> Also, the cell vectors-related procedure does not work with the ASE MD
> functions.
>
>
> GRADIENT NORM = 3903.37562 =
> 400.47787 PER ATOM
> IONIZATION POTENTIAL = 5.403926 EV
> HOMO LUMO ENERGIES (EV) = -5.404 -1.421
> NO. OF FILLED LEVELS = 94
> MOLECULAR WEIGHT = 436.8466 POINT GROUP:
> C1
>
> VOLUME OF UNIT CELL = 1000.000 CUBIC ANGSTROMS
>
> DENSITY = 0.725 GRAMS/CC
> A = 10.000 ANGSTROMS
> B = 10.000 ANGSTROMS
> C = 10.000 ANGSTROMS
> ALPHA = 90.000 DEGREES
> BETA = 90.000 DEGREES
> GAMMA = 90.000 DEGREES
>
> Pressure required to constrain translation vectors
> Tv( 96) Pressure: 16.32 GPa
> Tv( 97) Pressure: 52.10 GPa
> Tv( 98) Pressure: 136.90 GPa
>
>
>
> * Stress tensor in GPa using Voigt notation (xx, yy, zz, yz, xz,
> xy): -26.397 -22.259 -23.777 5.216 -4.175 -1.922*
>
>
> On Mon, Oct 28, 2024 at 10:43 AM Adam Jackson <a.j.jackson at physics.org>
> wrote:
>
>> No, the temperature is as I've set it. I guess the Mopac calculator simply
>> does not have its "get_stress()" function.
>>
>> This is correct, the ASE MOPAC calculator does not read stress from the
>> output file. Perhaps it was not present (or difficult to read reliably)
>> when it was written; I’ve just run a quick calculation on diamond and in
>> that case MOPAC 22 did write a stress tensor in Voigt form that we might be
>> able to work with. Perhaps somebody can create a Pull Request including a
>> unit test; we do have MOPAC 22 on the test runners so once set up it can be
>> checked frequently.
>>
>>
>>
>> I also found it necessary to set ‘GEOM-OK’ in the MOPAC parameters or the
>> calculation tends to complain about the unit cell and quit early.
>>
>>
>>
>> All the best,
>>
>> Adam
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20241028/cd80310d/attachment-0001.htm>
More information about the ase-users
mailing list