-
Notifications
You must be signed in to change notification settings - Fork 15
/
activation_relaxation_technique.py
181 lines (153 loc) · 6.55 KB
/
activation_relaxation_technique.py
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
# coding: utf-8
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
# Distributed under the terms of "New BSD License", see the LICENSE file.
import numpy as np
from pyiron_atomistics.atomistics.job.interactivewrapper import (
InteractiveWrapper,
ReferenceJobOutput,
)
from pyiron_base import DataContainer
__author__ = "Osamu Waseda"
__copyright__ = (
"Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH "
"- Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Osamu Waseda"
__email__ = "waseda@mpie.de"
__status__ = "development"
__date__ = "Sep 1, 2018"
class ARTInteractive(object):
def __init__(self, art_id, direction, gamma=0.1, fix_layer=False, non_art_id=None):
if int(art_id) != art_id or art_id < 0:
raise ValueError("art_id must be a posive integer")
if len(direction) != 3 or np.isclose(np.linalg.norm(direction), 0):
raise ValueError("direction must be a finite 3d vector")
if gamma < 0:
raise ValueError("gamma must be a positive float")
if fix_layer and non_art_id is not None:
raise ValueError("fix_layer and non_art_id cannot be set at the same time")
self.art_id = art_id
self.direction = direction
self.gamma = gamma
self.non_art_id = non_art_id
if non_art_id is not None:
self.non_art_id = np.array([non_art_id]).flatten()
self.fix_layer = fix_layer
@property
def _R(self):
value = np.array(self.direction)
value = value / np.linalg.norm(value)
return np.outer(value, value)
def get_forces(self, f_in):
f = np.array(f_in)
if len(f.shape) == 2:
f = np.array([f])
if self.non_art_id is None:
self.non_art_id = np.arange(len(f[0])) != self.art_id
f_art = (1.0 + self.gamma) * np.einsum("ij,nj->ni", self._R, f[:, self.art_id])
if self.fix_layer:
f[:, self.non_art_id] = np.einsum(
"nmj,ij->nmi", f[:, self.non_art_id], np.identity(3) - self._R
)
else:
f[:, self.non_art_id] += f_art[:, np.newaxis, :] / np.sum(
self.non_art_id != False
)
f[:, self.art_id] -= f_art
return f.reshape(np.array(f_in).shape)
class ART(InteractiveWrapper):
"""
Apply an artificial force according to the Activation Relaxation Technique (ART)
DOI:https://doi.org/10.1103/PhysRevE.57.2419
The applied force f_art is calculated from the original force f by:
f_art = f-(1+gamma)*np.dot(n,f)*n
where gamma is a parameter to be determined in the input (default: 0.1) and n is
the direction along which the force is reversed (3d-vector). In order to homogenize
the total force in the system, f_art is distributed among atoms specified by
non_art_id (default: all atoms), or if fix_layer is defined, the forces acting on
all the other atoms along the direction n are cancelled. Note: Since the energy
is not compatible with the forces anymore, structure optimization methods which
rely on the energy variation (such as conjugate gradient) cannot/should not be used.
Input:
- art_id (int): atom id on which ART force is applied
- direction (list/numpy.ndarray): direction along which force is reversed
- gamma (float): prefactor for force inversion. v.s.
- non_art_id (list/numpy.ndarray): list of atoms to be used for the force cancellation
- fix_layer (bool): whether or not to fix all other atoms on the layer perpendicular
to the direction along which force is reversed.
Example:
# Structure creation
>>> vacancy_position = structure.positions[0]
>>> del structure[0]
>>> neighbors = structure.get_neighborhood(vacancy_position)
>>> direction = neighbors.vecs[0]
>>> art_id = neighbors.indices[0]
>>> structure.positions[art_id] -= 0.5*direction
# Job creation
>>> some_atomistic_job.structure = structure
>>> art = ART(job_name='art')
>>> art.ref_job = some_atomistic_job
>>> art.input.art_id = art_id
>>> art.input.direction = direction
>>> some_minimizer.ref_job = art
>>> some_minimizer.run()
"""
def __init__(self, project, job_name):
super(ART, self).__init__(project, job_name)
self.input = DataContainer(table_name="custom_dict")
self.input.gamma = 0.1
self.input.fix_layer = False
self.input.non_art_id = None
self.input.art_id = None
self.input.direction = None
self.output = ARTIntOutput(job=self)
self.server.run_mode.interactive = True
self._interactive_interface = None
self._art = None
def set_input_to_read_only(self):
"""
This function enforces read-only mode for the input classes, but it has to be implement in the individual
classes.
"""
self.input.read_only = True
@property
def art(self):
if self._art is None:
self._art = ARTInteractive(
art_id=self.input.art_id,
direction=self.input.direction,
gamma=self.input.gamma,
fix_layer=self.input.fix_layer,
non_art_id=self.input.non_art_id,
)
return self._art
def run_if_interactive(self):
self._logger.debug("art status: " + str(self.status))
if not self.status.running:
self.ref_job_initialize()
self.status.running = True
if self.ref_job.server.run_mode.interactive:
self.ref_job.run()
else:
self.ref_job.run(run_again=True)
self._logger.debug("art status: " + str(self.status))
def interactive_forces_getter(self):
return self.art.get_forces(self.ref_job.output.forces[-1])
def validate_ready_to_run(self):
"""
check whether art_id and direction are set
"""
if self.input.art_id is None or self.input.direction is None:
raise AssertionError("art_id and/or direction not set")
def interactive_close(self):
self.status.collect = True
if self.ref_job.server.run_mode.interactive:
self.ref_job.interactive_close()
self.status.finished = True
class ARTIntOutput(ReferenceJobOutput):
def __init__(self, job):
super(ARTIntOutput, self).__init__(job=job)
@property
def forces(self):
return self._job.art.get_forces(self._job.ref_job.output.forces)