File: ga_bulk_relax.py

package info (click to toggle)
python-ase 3.26.0-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 15,484 kB
  • sloc: python: 148,112; xml: 2,728; makefile: 110; javascript: 47
file content (58 lines) | stat: -rw-r--r-- 1,870 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
from ase.build import niggli_reduce
from ase.calculators.singlepoint import SinglePointCalculator
from ase.constraints import ExpCellFilter
from ase.ga import set_raw_score
from ase.optimize import FIRE

try:
    from asap3 import EMT
except ImportError:
    from ase.calculators.emt import EMT


def finalize(atoms, energy=None, forces=None, stress=None):
    # Finalizes the atoms by attaching a SinglePointCalculator
    # and setting the raw score as the negative of the total energy
    atoms.wrap()
    calc = SinglePointCalculator(
        atoms, energy=energy, forces=forces, stress=stress
    )
    atoms.calc = calc
    raw_score = -atoms.get_potential_energy()
    set_raw_score(atoms, raw_score)


def relax(atoms, cellbounds=None):
    # Performs a variable-cell relaxation of the structure
    calc = EMT()
    atoms.calc = calc

    converged = False
    niter = 0
    while not converged and niter < 10:
        if cellbounds is not None:
            cell = atoms.get_cell()
            if not cellbounds.is_within_bounds(cell):
                niggli_reduce(atoms)
            cell = atoms.get_cell()
            if not cellbounds.is_within_bounds(cell):
                # Niggli reduction did not bring the unit cell
                # within the specified bounds; this candidate should
                # be discarded so we set an absurdly high energy
                finalize(atoms, 1e9)
                return

        ecf = ExpCellFilter(atoms)
        dyn = FIRE(ecf, maxmove=0.2, logfile=None, trajectory=None)
        dyn.run(fmax=1e-3, steps=100)

        converged = dyn.converged()
        niter += 1

    dyn = FIRE(atoms, maxmove=0.2, logfile=None, trajectory=None)
    dyn.run(fmax=1e-2, steps=100)

    e = atoms.get_potential_energy()
    f = atoms.get_forces()
    s = atoms.get_stress()
    finalize(atoms, energy=e, forces=f, stress=s)