diff --git a/model_zoo/research/hpc/sponge/main.py b/model_zoo/research/hpc/sponge/main.py index 9f37635f6c8..503946d8370 100644 --- a/model_zoo/research/hpc/sponge/main.py +++ b/model_zoo/research/hpc/sponge/main.py @@ -16,14 +16,14 @@ import argparse import time -from src.simulation import Simulation -from src.mdnn import Mdnn, TransCrdToCV import mindspore.context as context from mindspore import Tensor from mindspore import load_checkpoint +from src.mdnn import Mdnn, TransCrdToCV +from src.simulation import Simulation parser = argparse.ArgumentParser(description='SPONGE Controller') -parser.add_argument('--i', type=str, default=None, help='Input file') +parser.add_argument('--i', type=str, default=None, help='Input .in file') parser.add_argument('--amber_parm', type=str, default=None, help='Paramter file in AMBER type') parser.add_argument('--c', type=str, default=None, help='Initial coordinates file') parser.add_argument('--r', type=str, default="restrt", help='') @@ -36,6 +36,7 @@ parser.add_argument('--checkpoint', type=str, default="", help='Checkpoint file' args_opt = parser.parse_args() context.set_context(mode=context.GRAPH_MODE, device_target="GPU", device_id=args_opt.device_id, save_graphs=False) +# context.set_context(mode=context.PYNATIVE_MODE, device_target="GPU", device_id=args_opt.device_id, save_graphs=False) if __name__ == "__main__": simulation = Simulation(args_opt) @@ -53,7 +54,8 @@ if __name__ == "__main__": if steps == simulation.md_info.step_limit - 1: print_step = 0 temperature, total_potential_energy, sigma_of_bond_ene, sigma_of_angle_ene, sigma_of_dihedral_ene, \ - nb14_lj_energy_sum, nb14_cf_energy_sum, LJ_energy_sum, ee_ene, _ = simulation(Tensor(steps), Tensor(print_step)) + nb14_lj_energy_sum, nb14_cf_energy_sum, LJ_energy_sum, ee_ene, _, _, _, _ = simulation(Tensor(steps), + Tensor(print_step)) if steps == 0: compiler_time = time.time() diff --git a/model_zoo/research/hpc/sponge/src/angle.py b/model_zoo/research/hpc/sponge/src/angle.py index 38a1e4f3a79..a8e90dd4aae 100644 --- a/model_zoo/research/hpc/sponge/src/angle.py +++ b/model_zoo/research/hpc/sponge/src/angle.py @@ -13,12 +13,46 @@ # limitations under the License. # ============================================================================ '''Angle''' + + class Angle: '''Angle''' + def __init__(self, controller): + self.module_name = "angle" + self.h_atom_a = [] + self.h_atom_b = [] + self.h_atom_c = [] + self.h_angle_k = [] + self.h_angle_theta0 = [] + self.angle_numbers = 0 if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) + self.is_initialized = 1 + else: + self.read_in_file(controller) + + def read_in_file(self, controller): + """read_in_file""" + print("START INITIALIZING ANGLE:") + name = self.module_name + "_in_file" + if name in controller.Command_Set: + path = controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.angle_numbers = int(context[0].strip()) + print(" angle_numbers is ", self.angle_numbers) + for i in range(self.angle_numbers): + val = list(map(float, context[i + 1].strip().split())) + self.h_atom_a.append(int(val[0])) + self.h_atom_b.append(int(val[1])) + self.h_atom_c.append(int(val[2])) + self.h_angle_k.append(val[3]) + self.h_angle_theta0.append(val[4]) + self.is_initialized = 1 + file.close() + print("END INITIALIZING ANGLE") def read_information_from_amberfile(self, file_path): '''read amber file''' @@ -64,9 +98,9 @@ class Angle: information.extend(value) count += len(value) for _ in range(self.angle_with_H_numbers): - self.h_atom_a[angle_count] = information[angle_count * 4 + 0] / 3 - self.h_atom_b[angle_count] = information[angle_count * 4 + 1] / 3 - self.h_atom_c[angle_count] = information[angle_count * 4 + 2] / 3 + self.h_atom_a[angle_count] = int(information[angle_count * 4 + 0] / 3) + self.h_atom_b[angle_count] = int(information[angle_count * 4 + 1] / 3) + self.h_atom_c[angle_count] = int(information[angle_count * 4 + 2] / 3) self.h_type[angle_count] = information[angle_count * 4 + 3] - 1 angle_count += 1 @@ -86,9 +120,9 @@ class Angle: information.extend(value) count += len(value) for _ in range(self.angle_without_H_numbers): - self.h_atom_a[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 0] / 3 - self.h_atom_b[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 1] / 3 - self.h_atom_c[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 2] / 3 + self.h_atom_a[angle_count] = int(information[(angle_count - self.angle_with_H_numbers) * 4 + 0] / 3) + self.h_atom_b[angle_count] = int(information[(angle_count - self.angle_with_H_numbers) * 4 + 1] / 3) + self.h_atom_c[angle_count] = int(information[(angle_count - self.angle_with_H_numbers) * 4 + 2] / 3) self.h_type[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 3] - 1 angle_count += 1 break diff --git a/model_zoo/research/hpc/sponge/src/bd_baro.py b/model_zoo/research/hpc/sponge/src/bd_baro.py new file mode 100644 index 00000000000..b8a13780fb7 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/bd_baro.py @@ -0,0 +1,47 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''Angle''' + + +class BD_BARO: + '''Angle''' + + def __init__(self, controller, target_pressure, box_length, mode): + self.constant_pres_convertion = 6.946827162543585e4 + print("START INITIALIZING MC BAROSTAT:\n") + self.module_name = "bd_baro" + self.is_controller_printf_initialized = 0 + print(" The target pressure is %.2f bar\n" % (target_pressure * self.constant_pres_convertion)) + self.V0 = box_length[0] * box_length[1] * box_length[2] + self.newV = self.V0 + self.dt = 1e-3 if "dt" not in controller.Command_Set else float(controller.Command_Set["dt"]) + print(" The dt is %f ps\n" % (self.dt)) + self.taup = 1.0 if "tau" not in controller.Command_Set else float(controller.Command_Set["tau"]) + print(" The time constant tau is %f ps\n" % (self.taup)) + self.compressibility = 4.5e-5 if "compressibility" not in controller.Command_Set else float( + controller.Command_Set["compressibility"]) + print(" The compressibility constant is %f bar^-1\n" % (self.compressibility)) + self.update_interval = 10 if "update_interval" not in controller.Command_Set else int( + controller.Command_Set["update_interval"]) + print(" The update_interval is %d\n" % (self.update_interval)) + self.system_reinitializing_count = 0 + if mode == 2 and controller.Command_Set["barostat"] == "berendsen": + self.is_initialized = 1 + else: + self.is_initialized = 0 + + if self.is_initialized and not self.is_controller_printf_initialized: + self.is_controller_printf_initialized = 1 + print("END INITIALIZING BERENDSEN BAROSTATn") diff --git a/model_zoo/research/hpc/sponge/src/bond.py b/model_zoo/research/hpc/sponge/src/bond.py index 4cc5b659bd4..e0287f115e5 100644 --- a/model_zoo/research/hpc/sponge/src/bond.py +++ b/model_zoo/research/hpc/sponge/src/bond.py @@ -13,15 +13,45 @@ # limitations under the License. # ============================================================================ '''Bond''' + + class Bond: '''Bond''' - def __init__(self, controller, md_info): - - self.atom_numbers = md_info.atom_numbers + def __init__(self, controller): + self.module_name = "bond" + self.h_atom_a = [] + self.h_atom_b = [] + self.h_k = [] + self.h_r0 = [] + self.bond_numbers = 0 + self.is_initialized = 0 if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) + self.is_initialized = 1 + else: + self.read_in_file(controller) + + def read_in_file(self, controller): + """read_in_file""" + print("START INITIALIZING BOND:") + name = self.module_name + "_in_file" + if name in controller.Command_Set: + path = controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.bond_numbers = int(context[0].strip()) + print(" bond_numbers is ", self.bond_numbers) + for i in range(self.bond_numbers): + val = list(map(float, context[i + 1].strip().split())) + self.h_atom_a.append(int(val[0])) + self.h_atom_b.append(int(val[1])) + self.h_k.append(val[2]) + self.h_r0.append(val[3]) + self.is_initialized = 1 + file.close() + print("END INITIALIZING BOND") def read_information_from_amberfile(self, file_path): '''read amber file''' @@ -103,8 +133,8 @@ class Bond: count += len(value) for i in range(self.bond_with_hydrogen): - self.h_atom_a[i] = information[3 * i + 0] / 3 - self.h_atom_b[i] = information[3 * i + 1] / 3 + self.h_atom_a[i] = int(information[3 * i + 0] / 3) + self.h_atom_b[i] = int(information[3 * i + 1] / 3) tmpi = information[3 * i + 2] - 1 self.h_k[i] = self.bond_type_k[tmpi] self.h_r0[i] = self.bond_type_r[tmpi] @@ -126,8 +156,8 @@ class Bond: count += len(value) for i in range(self.bond_with_hydrogen, self.bond_numbers): - self.h_atom_a[i] = information[3 * (i - self.bond_with_hydrogen) + 0] / 3 - self.h_atom_b[i] = information[3 * (i - self.bond_with_hydrogen) + 1] / 3 + self.h_atom_a[i] = int(information[3 * (i - self.bond_with_hydrogen) + 0] / 3) + self.h_atom_b[i] = int(information[3 * (i - self.bond_with_hydrogen) + 1] / 3) tmpi = information[3 * (i - self.bond_with_hydrogen) + 2] - 1 self.h_k[i] = self.bond_type_k[tmpi] self.h_r0[i] = self.bond_type_r[tmpi] diff --git a/model_zoo/research/hpc/sponge/src/crd_molecular_map.py b/model_zoo/research/hpc/sponge/src/crd_molecular_map.py new file mode 100644 index 00000000000..dfc95c12261 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/crd_molecular_map.py @@ -0,0 +1,90 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''Coordinate Molecular Map''' +import collections +import numpy as np + + +class CoordinateMolecularMap: + '''Coordinate Molecular Map''' + + def __init__(self, atom_numbers, box_length, crd, exclude_numbers, exclude_length, exclude_start, exclude_list): + self.module_name = "crd_mole_wrap" + self.atom_numbers = atom_numbers + self.box_length = box_length + self.crd = crd + # self.coordinate = crd + self.exclude_numbers = exclude_numbers + self.exclude_length = exclude_length + self.exclude_start = exclude_start + self.exclude_list = exclude_list + print("START INITIALIZING Coordinate Molecular Map:\n") + self.h_box_map_times = np.zeros([self.atom_numbers, 3]) + self.h_old_crd = self.crd + self.h_nowrap_crd = self.crd + self.move_crd_nearest_from_exclusions_host() + self.is_initialized = 1 + print("END INITIALIZING Coordinate Molecular Map\n") + + def move_crd_nearest_from_exclusions_host(self): + '''move_crd_nearest_from_exclusions_host''' + edge_numbers = 2 * self.exclude_numbers + visited = [0] * self.atom_numbers + first_edge = [-1] * self.atom_numbers + edges = [0] * edge_numbers + edge_next = [0] * edge_numbers + atom_i, atom_j, edge_count = 0, 0, 0 + for i in range(self.atom_numbers): + atom_i = i + for j in range(self.exclude_start[i] + self.exclude_length[i] - 1, self.exclude_start[i] - 1, -1): + atom_j = self.exclude_list[j] + edge_next[edge_count] = first_edge[atom_i] + first_edge[atom_i] = edge_count + edges[edge_count] = atom_j + edge_count += 1 + edge_next[edge_count] = first_edge[atom_j] + first_edge[atom_j] = edge_count + edges[edge_count] = atom_i + edge_count += 1 + queue = collections.deque() + + for i in range(self.atom_numbers): + if not visited[i]: + visited[i] = 1 + queue.append(i) + atom_front = i + while queue: + atom = queue[0] + queue.popleft() + self.h_box_map_times[atom][0] = int( + (self.crd[atom_front][0] - self.crd[atom][0]) / self.box_length[0] \ + + 0.5) + self.h_box_map_times[atom][1] = int( + (self.crd[atom_front][1] - self.crd[atom][1]) / self.box_length[1] \ + + 0.5) + self.h_box_map_times[atom][2] = int( + (self.crd[atom_front][2] - self.crd[atom][2]) / self.box_length[1] \ + + 0.5) + self.crd[atom][0] = self.crd[atom][0] + self.h_box_map_times[atom][0] * self.box_length[0] + self.crd[atom][1] = self.crd[atom][1] + self.h_box_map_times[atom][1] * self.box_length[1] + self.crd[atom][2] = self.crd[atom][2] + self.h_box_map_times[atom][2] * self.box_length[2] + edge_count = first_edge[atom] + atom_front = atom + while edge_count is not -1: + atom = edges[edge_count] + if not visited[atom]: + queue.append(atom) + visited[atom] = 1 + edge_count = edge_next[edge_count] diff --git a/model_zoo/research/hpc/sponge/src/dihedral.py b/model_zoo/research/hpc/sponge/src/dihedral.py index 2d06c0e3b13..0eed5f9a8a0 100644 --- a/model_zoo/research/hpc/sponge/src/dihedral.py +++ b/model_zoo/research/hpc/sponge/src/dihedral.py @@ -18,11 +18,52 @@ import math class Dihedral: '''Dihedral''' + def __init__(self, controller): self.CONSTANT_Pi = 3.1415926535897932 + self.module_name = "dihedral" + self.h_atom_a = [] + self.h_atom_b = [] + self.h_atom_c = [] + self.h_atom_d = [] + self.h_ipn = [] + self.h_pn = [] + self.h_pk = [] + self.h_gamc = [] + self.h_gams = [] + self.dihedral_numbers = 0 if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) + self.is_initialized = 1 + else: + self.read_in_file(controller) + + def read_in_file(self, controller): + """read_in_file""" + print("START INITIALIZING DIHEDRAL:") + name = self.module_name + "_in_file" + if name in controller.Command_Set: + path = controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.dihedral_numbers = int(context[0].strip()) + print(" dihedral_numbers is ", self.dihedral_numbers) + for i in range(self.dihedral_numbers): + val = list(map(float, context[i + 1].strip().split())) + self.h_atom_a.append(int(val[0])) + self.h_atom_b.append(int(val[1])) + self.h_atom_c.append(int(val[2])) + self.h_atom_d.append(int(val[3])) + self.h_ipn.append(val[4]) + self.h_pn.append(val[4]) + self.h_pk.append(val[5]) + self.h_gamc.append(math.cos(val[6]) * val[5]) + self.h_gams.append(math.sin(val[6]) * val[5]) + + self.is_initialized = 1 + file.close() + print("END INITIALIZING DIHEDRAL") def read_information_from_amberfile(self, file_path): '''read amber file''' @@ -108,11 +149,11 @@ class Dihedral: self.h_atom_b = [0] * self.dihedral_numbers self.h_atom_c = [0] * self.dihedral_numbers self.h_atom_d = [0] * self.dihedral_numbers - self.pk = [] - self.gamc = [] - self.gams = [] - self.pn = [] - self.ipn = [] + self.h_pk = [] + self.h_gamc = [] + self.h_gams = [] + self.h_pn = [] + self.h_ipn = [] for idx, val in enumerate(context): if "%FLAG DIHEDRALS_INC_HYDROGEN" in val: count = 0 @@ -132,20 +173,20 @@ class Dihedral: self.h_atom_c[i] = information[i * 5 + 2] / 3 self.h_atom_d[i] = abs(information[i * 5 + 3] / 3) tmpi = information[i * 5 + 4] - 1 - self.pk.append(self.pk_type[tmpi]) + self.h_pk.append(self.pk_type[tmpi]) tmpf = self.phase_type[tmpi] if abs(tmpf - self.CONSTANT_Pi) <= 0.001: tmpf = self.CONSTANT_Pi tmpf2 = math.cos(tmpf) if abs(tmpf2) < 1e-6: tmpf2 = 0 - self.gamc.append(tmpf2 * self.pk[i]) + self.h_gamc.append(tmpf2 * self.h_pk[i]) tmpf2 = math.sin(tmpf) if abs(tmpf2) < 1e-6: tmpf2 = 0 - self.gams.append(tmpf2 * self.pk[i]) - self.pn.append(abs(self.pn_type[tmpi])) - self.ipn.append(int(self.pn[i] + 0.001)) + self.h_gams.append(tmpf2 * self.h_pk[i]) + self.h_pn.append(abs(self.pn_type[tmpi])) + self.h_ipn.append(int(self.h_pn[i] + 0.001)) break for idx, val in enumerate(context): if "%FLAG DIHEDRALS_WITHOUT_HYDROGEN" in val: @@ -166,20 +207,20 @@ class Dihedral: self.h_atom_c[i] = information[(i - self.dihedral_with_hydrogen) * 5 + 2] / 3 self.h_atom_d[i] = abs(information[(i - self.dihedral_with_hydrogen) * 5 + 3] / 3) tmpi = information[(i - self.dihedral_with_hydrogen) * 5 + 4] - 1 - self.pk.append(self.pk_type[tmpi]) + self.h_pk.append(self.pk_type[tmpi]) tmpf = self.phase_type[tmpi] if abs(tmpf - self.CONSTANT_Pi) <= 0.001: tmpf = self.CONSTANT_Pi tmpf2 = math.cos(tmpf) if abs(tmpf2) < 1e-6: tmpf2 = 0 - self.gamc.append(tmpf2 * self.pk[i]) + self.h_gamc.append(tmpf2 * self.h_pk[i]) tmpf2 = math.sin(tmpf) if abs(tmpf2) < 1e-6: tmpf2 = 0 - self.gams.append(tmpf2 * self.pk[i]) - self.pn.append(abs(self.pn_type[tmpi])) - self.ipn.append(int(self.pn[i] + 0.001)) + self.h_gams.append(tmpf2 * self.h_pk[i]) + self.h_pn.append(abs(self.pn_type[tmpi])) + self.h_ipn.append(int(self.h_pn[i] + 0.001)) break for i in range(self.dihedral_numbers): if self.h_atom_c[i] < 0: diff --git a/model_zoo/research/hpc/sponge/src/langevin_liujian_md.py b/model_zoo/research/hpc/sponge/src/langevin_liujian_md.py index 0f25929f9d5..6552f2b23df 100644 --- a/model_zoo/research/hpc/sponge/src/langevin_liujian_md.py +++ b/model_zoo/research/hpc/sponge/src/langevin_liujian_md.py @@ -20,37 +20,72 @@ import numpy as np class Langevin_Liujian: '''LagevinLiuJian''' + def __init__(self, controller, atom_numbers): + self.module_name = "langevin_liu" self.atom_numbers = atom_numbers + self.h_mass = [] + print("START INITIALIZING LANGEVIN_LIU DYNAMICS:") if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) - + else: + self.read_mass_file(controller) self.CONSTANT_TIME_CONVERTION = 20.455 self.CONSTANT_kB = 0.00198716 self.target_temperature = 300.0 if "target_temperature" not in controller.Command_Set else float( controller.Command_Set["target_temperature"]) - self.gamma_ln = 1.0 if "langevin_gamma" not in controller.Command_Set else float( - controller.Command_Set["langevin_gamma"]) - self.rand_seed = 1 if "langevin_seed" not in controller.Command_Set else float( - controller.Command_Set["langevin_seed"]) - self.max_velocity = 10000.0 if "velocity_max" not in controller.Command_Set else float( - controller.Command_Set["velocity_max"]) - assert self.max_velocity > 0 - print("target temperature is ", self.target_temperature) - print("friction coefficient is ", self.gamma_ln, "ps^-1") - print("random seed is ", self.rand_seed) - self.dt = float(controller.Command_Set["dt"]) - self.dt *= self.CONSTANT_TIME_CONVERTION + self.gamma_ln = 1.0 + if "gamma" in controller.Command_Set: + self.gamma_ln = float(controller.Command_Set["gamma"]) + if "langevin_liu_gamma" in controller.Command_Set: + self.gamma_ln = float(controller.Command_Set["langevin_liu_gamma"]) + print(" langevin_liu_gamma is ", self.gamma_ln) + + self.random_seed = 1 if "seed" not in controller.Command_Set else int( + controller.Command_Set["seed"]) + + print(" target temperature is {} K".format(self.target_temperature)) + print(" friction coefficient is {} ps^-1".format(self.gamma_ln)) + print(" random seed is ", self.random_seed) + self.dt = 0.001 if "dt" not in controller.Command_Set else float( + controller.Command_Set["dt"]) * self.CONSTANT_TIME_CONVERTION self.half_dt = 0.5 * self.dt - self.rand_state = np.float32(np.zeros([math.ceil(3 * self.atom_numbers / 4.0) * 16,])) + + self.float4_numbers = math.ceil(3 * self.atom_numbers / 4.0) + self.rand_state = np.float32(np.zeros([self.float4_numbers * 16,])) self.gamma_ln = self.gamma_ln / self.CONSTANT_TIME_CONVERTION self.exp_gamma = math.exp(-1 * self.gamma_ln * self.dt) self.sqrt_gamma = math.sqrt((1. - self.exp_gamma * self.exp_gamma) * self.target_temperature * self.CONSTANT_kB) self.h_sqrt_mass = [0] * self.atom_numbers for i in range(self.atom_numbers): - self.h_sqrt_mass[i] = self.sqrt_gamma * math.sqrt(1. / self.h_mass[i]) + self.h_sqrt_mass[i] = self.sqrt_gamma * math.sqrt(1. / self.h_mass[i]) if self.h_mass[i] != 0 else 0 + + self.max_velocity = 0 + if "velocity_max" in controller.Command_Set: + self.max_velocity = float(controller.Command_Set["velocity_max"]) + if "langevin_liu_velocity_max" in controller.Command_Set: + self.max_velocity = float(controller.Command_Set["langevin_liu_velocity_max"]) + print(" max velocity is ", self.max_velocity) + + self.h_mass_inverse = [0] * self.atom_numbers + for i in range(self.atom_numbers): + self.h_mass_inverse[i] = 1. / self.h_mass[i] if self.h_mass[i] != 0 else 0 + + self.is_initialized = 1 + + print("END INITIALIZING LANGEVIN_LIU DYNAMICS") + + def read_mass_file(self, controller): + if "mass_in_file" in controller.Command_Set: + path = controller.Command_Set["mass_in_file"] + file = open(path, 'r') + context = file.readlines() + for idx, val in enumerate(context): + if idx > 0: + self.h_mass.append(float(val.strip())) + file.close() def read_information_from_amberfile(self, file_path): '''read amber file''' diff --git a/model_zoo/research/hpc/sponge/src/lennard_jones.py b/model_zoo/research/hpc/sponge/src/lennard_jones.py index b7617c11d66..4b92affb7c7 100644 --- a/model_zoo/research/hpc/sponge/src/lennard_jones.py +++ b/model_zoo/research/hpc/sponge/src/lennard_jones.py @@ -13,12 +13,95 @@ # limitations under the License. # ============================================================================ '''Lennard Jones''' +import mindspore.common.dtype as mstype +from mindspore import Tensor +from mindspore.ops import operations as P + + class Lennard_Jones_Information: '''Lennard Jones''' - def __init__(self, controller): + + def __init__(self, controller, cutoff, box_length): + self.module_name = "LJ" + self.is_initialized = 0 + self.CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + self.CONSTANT_Pi = 3.1415926535897932 + self.cutoff = cutoff + self.box_length = box_length + if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) + self.is_initialized = 1 + else: + self.read_in_file(controller) + + if self.is_initialized: + self.totalc6get = P.totalc6get(self.atom_numbers) + self.read_information() + + def read_in_file(self, controller): + """read_in_file""" + print("START INITIALIZING LENNADR JONES INFORMATION:") + name = self.module_name + "_in_file" + # print("read_in_file " + name) + if name in controller.Command_Set: + path = controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.atom_numbers, self.atom_type_numbers = map(int, context[0].strip().split()) + print(" atom_numbers is ", self.atom_numbers) + print(" atom_LJ_type_number is ", self.atom_type_numbers) + self.pair_type_numbers = self.atom_type_numbers * (self.atom_type_numbers + 1) / 2 + self.h_LJ_A = [] + self.h_LJ_B = [] + self.h_atom_LJ_type = [] + startidx = 1 + count = 0 + print(startidx) + while count < self.atom_type_numbers: + if context[startidx].strip(): + val = list(map(float, context[startidx].strip().split())) + # print(val) + count += 1 + self.h_LJ_A.extend(val) + startidx += 1 + assert len(self.h_LJ_A) == self.pair_type_numbers + self.h_LJ_A = [x * 12.0 for x in self.h_LJ_A] + + count = 0 + print(startidx) + while count < self.atom_type_numbers: + if context[startidx].strip(): + val = list(map(float, context[startidx].strip().split())) + # print(val) + count += 1 + self.h_LJ_B.extend(val) + startidx += 1 + assert len(self.h_LJ_B) == self.pair_type_numbers + self.h_LJ_B = [x * 6.0 for x in self.h_LJ_B] + for idx, val in enumerate(context): + if idx > startidx: + self.h_atom_LJ_type.append(int(val.strip())) + file.close() + self.is_initialized = 1 + print("END INITIALIZING LENNADR JONES INFORMATION") + + def read_information(self): + """read_information""" + self.uint_dr_to_dr_cof = [1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[0], + 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[1], + 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[2]] + print("copy lj type to new crd") + self.atom_LJ_type = Tensor(self.h_atom_LJ_type, mstype.int32) + self.LJ_B = Tensor(self.h_LJ_B, mstype.float32) + self.factor = self.totalc6get(self.atom_LJ_type, self.LJ_B) + print(" factor is: ", self.factor) + self.long_range_factor = float(self.factor.asnumpy()) + self.long_range_factor *= -2.0 / 3.0 * self.CONSTANT_Pi / self.cutoff / self.cutoff / self.cutoff / 6.0 + self.volume = self.box_length[0] * self.box_length[1] * self.box_length[1] + print(" long range correction factor is: ", self.long_range_factor) + print(" End initializing long range LJ correction") def read_information_from_amberfile(self, file_path): '''read amber file''' @@ -35,9 +118,9 @@ class Lennard_Jones_Information: self.atom_numbers = value[0] self.atom_type_numbers = value[1] self.pair_type_numbers = int( - self.atom_type_numbers * (self.atom_type_numbers + 1) / 2) # TODO 这个地方有问题啊 + self.atom_type_numbers * (self.atom_type_numbers + 1) / 2) # TODO break - self.atom_LJ_type = [0] * self.atom_numbers + self.h_atom_LJ_type = [0] * self.atom_numbers for idx, val in enumerate(context): if "%FLAG ATOM_TYPE_INDEX" in val: count = 0 @@ -52,9 +135,9 @@ class Lennard_Jones_Information: information.extend(value) count += len(value) for i in range(self.atom_numbers): - self.atom_LJ_type[i] = information[i] - 1 + self.h_atom_LJ_type[i] = information[i] - 1 break - self.LJ_A = [0] * self.pair_type_numbers + self.h_LJ_A = [0] * self.pair_type_numbers for idx, val in enumerate(context): if "%FLAG LENNARD_JONES_ACOEF" in val: count = 0 @@ -69,9 +152,9 @@ class Lennard_Jones_Information: information.extend(value) count += len(value) for i in range(self.pair_type_numbers): - self.LJ_A[i] = 12.0 * information[i] + self.h_LJ_A[i] = 12.0 * information[i] break - self.LJ_B = [0] * self.pair_type_numbers + self.h_LJ_B = [0] * self.pair_type_numbers for idx, val in enumerate(context): if "%FLAG LENNARD_JONES_BCOEF" in val: count = 0 @@ -86,5 +169,5 @@ class Lennard_Jones_Information: information.extend(value) count += len(value) for i in range(self.pair_type_numbers): - self.LJ_B[i] = 6.0 * information[i] + self.h_LJ_B[i] = 6.0 * information[i] break diff --git a/model_zoo/research/hpc/sponge/src/mc_baro.py b/model_zoo/research/hpc/sponge/src/mc_baro.py new file mode 100644 index 00000000000..6de764978e9 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/mc_baro.py @@ -0,0 +1,68 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''MC_BARO''' + + +class MC_BARO: + '''MC_BARO''' + + def __init__(self, controller, atom_numbers, target_pressure, boxlength, res_is_initialized, mode): + self.constant_pres_convertion = 6.946827162543585e4 + print("START INITIALIZING MC BAROSTAT:\n") + self.module_name = "mc_baro" + self.atom_numbers = atom_numbers + self.accept_rate = 0 + self.DeltaV_max = 200.0 + self.is_controller_printf_initialized = 0 + self.is_initialized = 0 + # initial + self.target_pressure = target_pressure + print(" The target pressure is %.2f bar\n" % (target_pressure * self.constant_pres_convertion)) + self.V0 = boxlength[0] * boxlength[1] * boxlength[2] + self.newV = self.V0 + self.mc_baro_initial_ratio = 0.01 if "initial_ratio" not in controller.Command_Set else float( + controller.Command_Set["initial_ratio"]) + self.DeltaV_max = self.mc_baro_initial_ratio * self.V0 + print(" The initial max volume to change is %f A^3\n" % (self.DeltaV_max)) + self.update_interval = 100 if "update_interval" not in controller.Command_Set else int( + controller.Command_Set["update_interval"]) + print(" The update_interval is %d\n" % (self.update_interval)) + self.check_interval = 20 if "check_interval" not in controller.Command_Set else int( + controller.Command_Set["check_interval"]) + print(" The check_interval is %d\n" % (self.check_interval)) + self.scale_coordinate_by_residue = res_is_initialized if "residue_scale" not in controller.Command_Set else int( + controller.Command_Set["residue_scale"]) + if self.scale_coordinate_by_residue == 1 and res_is_initialized == 0: + print( + " Warning: The residue is not initialized. Atom scale mode is set instead.\n") + self.scale_coordinate_by_residue = 0 + print(" The residue_scale is %d\n" % (self.scale_coordinate_by_residue)) + + self.system_reinitializing_count = 0 + self.accept_rate_low = 30.0 if "accept_rate_low" not in controller.Command_Set else int( + controller.Command_Set["accept_rate_low"]) + print(" The lowest accept rate is %.2f\n" % (self.accept_rate_low)) + self.accept_rate_high = 40.0 if "accept_rate_high" not in controller.Command_Set else int( + controller.Command_Set["accept_rate_high"]) + print(" The highest accept rate is %.2f\n" % (self.accept_rate_high)) + self.frc_backup = 0 # [atom_numbers, 3] + self.crd_backup = 0 # [atom_numbers, 3] + if mode == 2 and controller.Command_Set["barostat"] == "monte_carlo": + self.is_initialized = 1 + else: + self.is_initialized = 0 + if self.is_initialized and not self.is_controller_printf_initialized: + self.is_controller_printf_initialized = 1 + print("END INITIALIZING MC BAROSTAT\n") diff --git a/model_zoo/research/hpc/sponge/src/md_information.py b/model_zoo/research/hpc/sponge/src/md_information.py index f4dc2e26f17..263608b8e99 100644 --- a/model_zoo/research/hpc/sponge/src/md_information.py +++ b/model_zoo/research/hpc/sponge/src/md_information.py @@ -14,53 +14,206 @@ # ============================================================================ '''MD Information''' import numpy as np +from src.system_information import (periodic_box_condition_information, system_information, + non_bond_information, NVE_iteration, residue_information, trajectory_output) class md_information: '''MD Information''' + def __init__(self, controller): CONSTANT_TIME_CONVERTION = 20.455 - CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + self.md_task = controller.md_task - self.mode = 0 if "mode" not in controller.Command_Set else int(controller.Command_Set["mode"]) - self.dt = 0.001 * CONSTANT_TIME_CONVERTION if "dt" not in controller.Command_Set else float( - controller.Command_Set["dt"]) * CONSTANT_TIME_CONVERTION - self.skin = 2.0 if "skin" not in controller.Command_Set else float(controller.Command_Set["skin"]) - self.trans_vec = [self.skin, self.skin, self.skin] - self.trans_vec_minus = -1 * self.trans_vec - self.step_limit = 1000 if "step_limit" not in controller.Command_Set else int( - controller.Command_Set["step_limit"]) + self.netfrc = 0 if "net_force" not in controller.Command_Set else int(controller.Command_Set["net_force"]) self.ntwx = 1000 if "write_information_interval" not in controller.Command_Set else int( controller.Command_Set["write_information_interval"]) - self.ntce = self.step_limit + 1 if "calculate_energy_interval" not in controller.Command_Set else int( - controller.Command_Set["calculate_energy_interval"]) self.atom_numbers = 0 self.residue_numbers = 0 self.density = 0.0 self.lin_serial = [] self.h_res_start = [] self.h_res_end = [] + + self.h_charge = [] self.h_mass = [] self.h_mass_inverse = [] self.h_charge = [] + self.coordinate = [] + self.box_length = [] + self.vel = [] + self.crd = [] + self.velocity = [] + + self.mode = self.read_mode(controller) + # read dt + self.dt = 0.001 * CONSTANT_TIME_CONVERTION if "dt" not in controller.Command_Set else float( + controller.Command_Set["dt"]) * CONSTANT_TIME_CONVERTION + self.dt_in_ps = 0.001 if "dt" not in controller.Command_Set else float(controller.Command_Set["dt"]) if controller.amber_parm is not None: self.read_basic_system_information_from_amber_file(controller.amber_parm) - - if "amber_irest" in controller.Command_Set: - amber_irest = int(controller.Command_Set["amber_irest"]) if controller.initial_coordinates_file is not None: - self.read_basic_system_information_from_rst7(controller.initial_coordinates_file, amber_irest) + self.read_basic_system_information_from_rst7(controller.initial_coordinates_file) + else: + self.read_coordinate_and_velocity(controller) + self.read_mass(controller) + self.read_charge(controller) + self.crd = self.coordinate - self.crd_to_uint_crd_cof = [CONSTANT_UINT_MAX_FLOAT / self.box_length[0], - CONSTANT_UINT_MAX_FLOAT / self.box_length[1], - CONSTANT_UINT_MAX_FLOAT / self.box_length[2]] - self.uint_dr_to_dr_cof = [1.0 / self.crd_to_uint_crd_cof[0], 1.0 / self.crd_to_uint_crd_cof[1], - 1.0 / self.crd_to_uint_crd_cof[2]] - self.density *= 1e24 / 6.023e23 / (self.box_length[0] * self.box_length[1] * self.box_length[2]) + self.sys = system_information(controller, self) + self.nb = non_bond_information(controller, self) + self.output = trajectory_output(controller, self) + self.nve = NVE_iteration(controller, self) + self.res = residue_information(controller, self) + self.pbc = periodic_box_condition_information(controller, self.box_length) + + if not self.h_res_start: + self.h_res_start = self.res.h_res_start + self.h_res_end = self.res.h_res_end + self.residue_numbers = self.res.residue_numbers + + # Atom_Information_Initial + self.acc = np.zeros([self.atom_numbers, 3]) + self.frc = np.zeros([self.atom_numbers, 3]) + self.sys.freedom = 3 * self.atom_numbers + self.is_initialized = 1 self.velocity = np.reshape(np.asarray(self.velocity, np.float32), [self.atom_numbers, 3]) + self.step_limit = self.sys.step_limit + + def read_mode(self, controller): + """read_mode""" + if "mode" in controller.Command_Set: + if controller.Command_Set["mode"] in ["NVT", "nvt", "1"]: + print(" Mode set to NVT\n") + mode = 1 + elif controller.Command_Set["mode"] in ["NPT", "npt", "2"]: + print(" Mode set to NPT\n") + mode = 2 + elif controller.Command_Set["mode"] in ["Minimization", "minimization", "-1"]: + print(" Mode set to Energy Minimization\n") + mode = -1 + elif controller.Command_Set["mode"] in ["NVE", "nve", "0"]: + print(" Mode set to NVE\n") + mode = 0 + else: + print( + " Warning: Mode {} is not match. Set to NVE as default\n".format(controller.Command_Set["mode"])) + mode = 0 + else: + print(" Mode set to NVE as default\n") + mode = 0 + return mode + + def read_coordinate_in_file(self, path): + '''read coordinates file''' + file = open(path, 'r') + print(" Start reading coordinate_in_file:\n") + context = file.readlines() + atom_numbers = int(context[0].strip()) + if self.atom_numbers != 0: + if self.atom_numbers is not atom_numbers: + print(" Error: atom_numbers is not equal: ", atom_numbers, self.atom_numbers) + exit(1) + else: + self.atom_numbers = atom_numbers + print(" atom_numbers is ", self.atom_numbers) + + for idx in range(self.atom_numbers): + coord = list(map(float, context[idx + 1].strip().split())) + self.coordinate.append(coord) + + self.box_length = list(map(float, context[-1].strip().split()))[:3] + print(" box_length is: x: {}, y: {}, z: {}".format( + self.box_length[0], self.box_length[1], self.box_length[2])) + self.crd = self.coordinate + file.close() + + def read_velocity_in_file(self, path): + '''read velocity file''' + file = open(path, 'r') + print(" Start reading velocity_in_file:\n") + context = file.readlines() + for idx, val in enumerate(context): + if idx == 0: + atom_numbers = int(val.strip()) + if self.atom_numbers > 0 and atom_numbers != self.atom_numbers: + print(" Error: atom_numbers is not equal: %d %d\n", idx, self.atom_numbers) + exit(1) + else: + self.atom_numbers = atom_numbers + else: + vel = list(map(float, val.strip().split())) + self.velocity.append(vel) + self.vel = self.velocity + file.close() + + def read_coordinate_and_velocity(self, controller): + """read_coordinate_and_velocity""" + if "coordinate_in_file" in controller.Command_Set: + self.read_coordinate_in_file(controller.Command_Set["coordinate_in_file"]) + if "velocity_in_file" in controller.Command_Set: + self.read_velocity_in_file(controller.Command_Set["velocity_in_file"]) + else: + print(" Velocity is set to zero as default\n") + self.velocity = [0] * 3 * self.atom_numbers + + def read_mass(self, controller): + """read_mass""" + print(" Start reading mass:") + if "mass_in_file" in controller.Command_Set: + path = controller.Command_Set["mass_in_file"] + file = open(path, 'r') + self.total_mass = 0 + context = file.readlines() + for idx, val in enumerate(context): + if idx == 0: + atom_numbers = int(val.strip()) + if self.atom_numbers > 0 and (atom_numbers != self.atom_numbers): + print(" Error: atom_numbers is not equal: ", atom_numbers, self.atom_numbers) + exit(1) + else: + self.atom_numbers = atom_numbers + else: + mass = float(val.strip()) + self.h_mass.append(mass) + self.total_mass += mass + if mass == 0: + self.h_mass_inverse.append(0.0) + else: + self.h_mass_inverse.append(1 / mass) + file.close() + else: + print(" mass is set to 20 as default") + self.total_mass = 20 * self.atom_numbers + self.h_mass = [20] * self.atom_numbers + self.h_mass_inverse = [1 / 20] * self.atom_numbers + + print(" End reading mass") + + def read_charge(self, controller): + """read_charge""" + if "charge_in_file" in controller.Command_Set: + print(" Start reading charge:") + path = controller.Command_Set["charge_in_file"] + file = open(path, 'r') + context = file.readlines() + for idx, val in enumerate(context): + if idx == 0: + atom_numbers = int(val.strip()) + if self.atom_numbers > 0 and (atom_numbers != self.atom_numbers): + print(" Error: atom_numbers is not equal: %d %d\n", idx, self.atom_numbers) + exit(1) + else: + self.atom_numbers = atom_numbers + else: + self.h_charge.append(float(val.strip())) + file.close() + else: + self.h_charge = [0.0] * self.atom_numbers + print(" End reading charge") def read_basic_system_information_from_amber_file(self, path): '''read amber file''' @@ -137,11 +290,13 @@ class md_information: count += len(value) break - def read_basic_system_information_from_rst7(self, path, irest): + def read_basic_system_information_from_rst7(self, path): '''read rst7 file''' file = open(path, 'r') context = file.readlines() file.close() + x = context[1].strip().split() + irest = 1 if len(x) > 1 else 0 atom_numbers = int(context[1].strip().split()[0]) if atom_numbers != self.atom_numbers: print("ERROR") @@ -151,7 +306,7 @@ class md_information: count = 0 start_idx = 1 if irest == 1: - self.simulation_start_time = float(context[1].strip().split()[1]) + self.simulation_start_time = float(x[1]) while count <= 6 * self.atom_numbers + 3: start_idx += 1 value = list(map(float, context[start_idx].strip().split())) @@ -169,4 +324,6 @@ class md_information: self.coordinate = information[: 3 * self.atom_numbers] self.velocity = [0.0] * (3 * self.atom_numbers) self.box_length = information[3 * self.atom_numbers:3 * self.atom_numbers + 3] + self.coordinate = np.array(self.coordinate).reshape([-1, 3]) + self.velocity = np.array(self.velocity).reshape([-1, 3]) print("system size is ", self.box_length[0], self.box_length[1], self.box_length[2]) diff --git a/model_zoo/research/hpc/sponge/src/nb14.py b/model_zoo/research/hpc/sponge/src/nb14.py index 9c37ec79e02..b28f13645d8 100644 --- a/model_zoo/research/hpc/sponge/src/nb14.py +++ b/model_zoo/research/hpc/sponge/src/nb14.py @@ -13,21 +13,51 @@ # limitations under the License. # ============================================================================ '''NON BOND''' + + class NON_BOND_14: '''NON BOND''' - def __init__(self, controller, dihedral, atom_numbers): - self.dihedral_with_hydrogen = dihedral.dihedral_with_hydrogen - self.dihedral_numbers = dihedral.dihedral_numbers - self.dihedral_type_numbers = dihedral.dihedral_type_numbers - self.atom_numbers = atom_numbers + def __init__(self, controller, dihedral, atom_numbers): + self.module_name = "nb14" + self.atom_numbers = atom_numbers + self.h_atom_a = [] + self.h_atom_b = [] + self.h_lj_scale_factor = [] + self.h_cf_scale_factor = [] + self.nb14_numbers = 0 + self.is_initialized = 0 if controller.amber_parm is not None: + self.dihedral_with_hydrogen = dihedral.dihedral_with_hydrogen + self.dihedral_numbers = dihedral.dihedral_numbers + self.dihedral_type_numbers = dihedral.dihedral_type_numbers file_path = controller.amber_parm self.read_information_from_amberfile(file_path) - self.h_atom_a = self.h_atom_a[:self.nb14_numbers] - self.h_atom_b = self.h_atom_b[:self.nb14_numbers] - self.h_lj_scale_factor = self.h_lj_scale_factor[:self.nb14_numbers] - self.h_cf_scale_factor = self.h_cf_scale_factor[:self.nb14_numbers] + self.h_atom_a = self.h_atom_a[:self.nb14_numbers] + self.h_atom_b = self.h_atom_b[:self.nb14_numbers] + self.h_lj_scale_factor = self.h_lj_scale_factor[:self.nb14_numbers] + self.h_cf_scale_factor = self.h_cf_scale_factor[:self.nb14_numbers] + self.is_initialized = 1 + else: + self.read_in_file(controller) + + def read_in_file(self, controller): + """read_in_file""" + name = self.module_name + "_in_file" + if name in controller.Command_Set: + path = controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.nb14_numbers = int(context[0].strip()) + print(" non-bond 14 numbers is", self.nb14_numbers) + for i in range(self.nb14_numbers): + val = list(map(float, context[i + 1].strip().split())) + self.h_atom_a.append(int(val[0])) + self.h_atom_b.append(int(val[1])) + self.h_lj_scale_factor.append(val[2]) + self.h_cf_scale_factor.append(val[3]) + self.is_initialized = 1 + file.close() def read_information_from_amberfile(self, file_path): '''read amber file''' diff --git a/model_zoo/research/hpc/sponge/src/neighbor_list.py b/model_zoo/research/hpc/sponge/src/neighbor_list.py index 607f6d258c2..81c5868bd56 100644 --- a/model_zoo/research/hpc/sponge/src/neighbor_list.py +++ b/model_zoo/research/hpc/sponge/src/neighbor_list.py @@ -13,17 +13,24 @@ # limitations under the License. # ============================================================================ '''Neighbor List''' + + class neighbor_list: '''Neighbor List''' + def __init__(self, controller, atom_numbers, box_length): - self.refresh_interval = 20 if "neighbor_list_refresh_interval" not in controller.Command_Set else int( - controller.Command_Set["neighbor_list_refresh_interval"]) + self.CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + print("START INITIALIZING NEIGHBOR LIST:") + self.module_name = "neighbor_list" + self.refresh_interval = 20 if "refresh_interval" not in controller.Command_Set else int( + controller.Command_Set["refresh_interval"]) self.max_atom_in_grid_numbers = 64 if "max_atom_in_grid_numbers" not in controller.Command_Set else int( controller.Command_Set["max_atom_in_grid_numbers"]) self.max_neighbor_numbers = 800 if "max_neighbor_numbers" not in controller.Command_Set else int( controller.Command_Set["max_neighbor_numbers"]) + self.skin = 2.0 if "skin" not in controller.Command_Set else float(controller.Command_Set["skin"]) - self.cutoff = 10.0 if "cut" not in controller.Command_Set else float(controller.Command_Set["cut"]) + self.cutoff = 10.0 if "cutoff" not in controller.Command_Set else float(controller.Command_Set["cutoff"]) self.cutoff_square = self.cutoff * self.cutoff self.cutoff_with_skin = self.cutoff + self.skin self.half_cutoff_with_skin = 0.5 * self.cutoff_with_skin @@ -31,15 +38,17 @@ class neighbor_list: self.half_skin_square = 0.25 * self.skin * self.skin self.atom_numbers = atom_numbers self.box_length = box_length + self.update_volume() + + self.initial_neighbor_grid() + self.not_first_time = 0 + self.is_initialized = 1 + self.refresh_count = [0] if controller.amber_parm is not None: file_path = controller.amber_parm self.read_information_from_amberfile(file_path) - self.Initial_Neighbor_Grid() - self.not_first_time = 0 - self.refresh_count = [0] - def read_information_from_amberfile(self, file_path): '''read amber file''' file = open(file_path, 'r') @@ -117,20 +126,23 @@ class neighbor_list: self.excluded_list.extend(tmp_list) break - def Initial_Neighbor_Grid(self): + def initial_neighbor_grid(self): '''init neighbor grid''' half_cutoff = self.half_cutoff_with_skin self.Nx = int(self.box_length[0] / half_cutoff) self.Ny = int(self.box_length[1] / half_cutoff) self.Nz = int(self.box_length[2] / half_cutoff) self.grid_N = [self.Nx, self.Ny, self.Nz] - self.grid_length = [self.box_length[0] / self.Nx, self.box_length[1] / self.Ny, self.box_length[2] / self.Nz] + self.grid_length = [self.box_length[0] / self.Nx, + self.box_length[1] / self.Ny, + self.box_length[2] / self.Nz] self.grid_length_inverse = [1.0 / self.grid_length[0], 1.0 / self.grid_length[1], 1.0 / self.grid_length[2]] + self.Nxy = self.Nx * self.Ny self.grid_numbers = self.Nz * self.Nxy - self.atom_numbers_in_grid_bucket = [0] * self.grid_numbers self.bucket = [-1] * (self.grid_numbers * self.max_atom_in_grid_numbers) + self.pointer = [] temp_grid_serial = [0] * 125 for i in range(self.grid_numbers): @@ -160,3 +172,11 @@ class neighbor_list: count += 1 temp_grid_serial = sorted(temp_grid_serial) self.pointer.extend(temp_grid_serial) + + def update_volume(self): + self.quarter_crd_to_uint_crd_cof = [0.25 * self.CONSTANT_UINT_MAX_FLOAT / self.box_length[0], + 0.25 * self.CONSTANT_UINT_MAX_FLOAT / self.box_length[1], + 0.25 * self.CONSTANT_UINT_MAX_FLOAT / self.box_length[2]] + self.uint_dr_to_dr_cof = [1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[0], + 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[1], + 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length[2]] diff --git a/model_zoo/research/hpc/sponge/src/particle_mesh_ewald.py b/model_zoo/research/hpc/sponge/src/particle_mesh_ewald.py index fd7f20f0104..4b22137d045 100644 --- a/model_zoo/research/hpc/sponge/src/particle_mesh_ewald.py +++ b/model_zoo/research/hpc/sponge/src/particle_mesh_ewald.py @@ -19,23 +19,40 @@ import math class Particle_Mesh_Ewald(): '''PME''' def __init__(self, controller, md_info): - self.cutoff = 10.0 if "cut" not in controller.Command_Set else float(controller.Command_Set["cut"]) - self.tolerance = 0.00001 if "PME_Direct_Tolerance" not in controller.Command_Set else float( - controller.Command_Set["PME_Direct_Tolerance"]) + self.module_name = "PME" + self.CONSTANT_Pi = 3.1415926535897932 + self.cutoff = 10.0 if "cutoff" not in controller.Command_Set else float(controller.Command_Set["cutoff"]) + self.tolerance = 0.00001 if "Direct_Tolerance" not in controller.Command_Set else float( + controller.Command_Set["Direct_Tolerance"]) self.fftx = -1 if "fftx" not in controller.Command_Set else int(controller.Command_Set["fftx"]) self.ffty = -1 if "ffty" not in controller.Command_Set else int(controller.Command_Set["ffty"]) self.fftz = -1 if "fftz" not in controller.Command_Set else int(controller.Command_Set["fftz"]) self.atom_numbers = md_info.atom_numbers self.box_length = md_info.box_length + self.volume = self.box_length[0] * self.box_length[1] * self.box_length[1] + if self.fftx < 0: self.fftx = self.Get_Fft_Patameter(self.box_length[0]) if self.ffty < 0: self.ffty = self.Get_Fft_Patameter(self.box_length[1]) if self.fftz < 0: self.fftz = self.Get_Fft_Patameter(self.box_length[2]) + print(" fftx: ", self.fftx) + print(" ffty: ", self.ffty) + print(" fftz: ", self.fftz) + print("pme cutoff", self.cutoff) + print("pme tolerance", self.tolerance) + self.PME_Nall = self.fftx * self.ffty * self.fftz + self.PME_Nin = self.ffty * self.fftz + self.PME_Nfft = self.fftx * self.ffty * (int(self.fftz / 2) + 1) + self.PME_inverse_box_vector = [self.fftx / self.box_length[0], + self.ffty / self.box_length[1], + self.fftz / self.box_length[2]] self.beta = self.Get_Beta(self.cutoff, self.tolerance) + self.neutralizing_factor = -0.5 * self.CONSTANT_Pi / (self.beta * self.beta * self.volume) + self.is_initialized = 1 def Get_Beta(self, cutoff, tolerance): '''GET BETA''' diff --git a/model_zoo/research/hpc/sponge/src/restrain.py b/model_zoo/research/hpc/sponge/src/restrain.py new file mode 100644 index 00000000000..2021b10b4d9 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/restrain.py @@ -0,0 +1,66 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''restrain_information''' + + +class Restrain_Information: + '''Angle''' + + def __init__(self, controller, atom_numbers, crd): + self.module_name = "restrain" + self.atom_numbers = atom_numbers + self.crd = crd + self.controller = controller + self.weight = 100.0 if "weight" not in controller.Command_Set else float( + controller.Command_Set["weight"]) + print(" %s_weight is %.0f\n" % (self.module_name, self.weight)) + self.is_initialized = 0 + name = self.module_name + "_in_file" + if name in controller.Command_Set: + print("START INITIALIZING RESTRAIN\n") + path = controller.Command_Set[name] + self.read_in_file(path) + self.read_crd_ref() + self.is_initialized = 1 + print("END INITIALIZING RESTRAIN\n") + + def read_in_file(self, path): + """read_in_file""" + file = open(path, 'r') + context = file.readlines() + self.restrain_numbers = 0 + h_lists = [] + for _, val in enumerate(context): + h_lists.append(float(val.strip())) + self.restrain_numbers += 1 + print(" restrain_numbers is %d\n", self.restrain_numbers) + file.close() + self.restrain_ene = [0] * self.restrain_numbers + self.sum_of_restrain_ene = [0] + + def read_crd_ref(self): + """read_crd_ref""" + self.crd_ref = [] + if "coordinate" not in self.controller.Command_Set: + print(" restrain reference coordinate copy from input coordinate") + self.crd_ref = self.crd + else: + print(" reading restrain reference coordinate file") + file = open(self.controller.Command_Set["coordinate"], 'r') + context = file.readlines() + atom_numbers = int(context[0].strip()) + print(" atom_numbers is %d", atom_numbers) + for i in range(atom_numbers): + self.crd_ref.append(list(map(float, context[i + 1].strip().split()))) diff --git a/model_zoo/research/hpc/sponge/src/simple_constrain.py b/model_zoo/research/hpc/sponge/src/simple_constrain.py new file mode 100644 index 00000000000..c40dee00235 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/simple_constrain.py @@ -0,0 +1,203 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''Simple_Constarin''''' +import math + + +class Bond_Information: + def __init__(self): + self.bond_numbers = 0 + self.atom_a = [] + self.atom_b = [] + self.bond_r = [] + + +class Information: + '''information''' + + def __init__(self): + self.atom_numbers = 0 + self.dt = 0.001 + self.uint_dr_to_dr_cof = [] + self.quarter_crd_to_uint_crd_cof = [] + self.volume = 0.0 + self.exp_gamma = 1.0 + self.step_length = 1.0 + self.iteration_numbers = 25 + self.constrain_mass = 3 + + +class Simple_Constarin: + '''Simple_Constarin''' + + def __init__(self, controller, md_info, bond, angle, liujian_info): + self.module_name = "simple_constrain" + self.CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + self.controller = controller + self.md_info = md_info + self.bond = bond + self.angle = angle + self.liujian_info = liujian_info + self.bond_info = Bond_Information() + self.info = Information() + + self.atom_numbers = md_info.atom_numbers + self.bond_numbers = bond.bond_numbers + self.angle_numbers = angle.angle_numbers + + self.is_initialized = 0 + self.constrain_frc = [] + self.pair_virial = [] + self.virial = [] + self.test_uint_crd = [] + self.last_pair_dr = [] + self.dt_inverse = 0 + self.half_exp_gamma_plus_half = 0 + self.bond_constrain_pair_numbers = 0 + self.angle_constrain_pair_numbers = 0 + self.constrain_pair_numbers = 0 + self.h_bond_pair = [] + self.h_constrain_pair = [] + if "constrain_mode" in self.controller.Command_Set and \ + self.controller.Command_Set["constrain_mode"] == "simple_constrain": + print("START INITIALIZING SIMPLE CONSTRAIN:") + self.constrain_pair_numbers = 0 + self.add_hbond_to_constrain_pair(self.bond_numbers, bond.h_atom_a, bond.h_atom_b, bond.h_r0, md_info.h_mass) + self.add_hangle_to_constrain_pair(self.angle_numbers, angle.h_atom_a, angle.h_atom_b, + angle.h_atom_c, angle.h_angle_theta0, md_info.h_mass) + if liujian_info.is_initialized: + self.initial_simple_constrain(md_info.atom_numbers, md_info.dt, md_info.sys.box_length, + liujian_info.exp_gamma, 0, md_info.h_mass, md_info.sys.freedom) + else: + self.initial_simple_constrain(md_info.atom_numbers, md_info.dt, md_info.sys.box_length, + 1.0, md_info.mode == -1, md_info.h_mass, md_info.sys.freedom) + self.is_initialized = 1 + print("END INITIALIZING SIMPLE CONSTRAIN\n") + + def add_hbond_to_constrain_pair(self, bond_numbers, atom_a, atom_b, bond_r, atom_mass): + """add_hbond_to_constrain_pair""" + self.info.constrain_mass = 3.0 + name = self.module_name + "_in_file" + if name in self.controller.Command_Set: + self.info.constrain_mass = 0 + if "mass" in self.controller.Command_Set: + self.info.constrain_mass = float(self.controller.Command_Set["mass"]) + self.h_bond_pair = [] + s = 0 + for i in range(bond_numbers): + mass_a = atom_mass[atom_a[i]] + mass_b = atom_mass[atom_b[i]] + if (0 < mass_a < self.info.constrain_mass) or (0 < mass_b < self.info.constrain_mass): + constrain_k = atom_mass[atom_a[i]] * atom_mass[atom_b[i]] / \ + (atom_mass[atom_a[i]] + atom_mass[atom_b[i]]) + self.h_bond_pair.append([atom_a[i], atom_b[i], bond_r[i], constrain_k]) + s += 1 + self.bond_constrain_pair_numbers = s + self.bond_info.bond_numbers = bond_numbers + self.bond_info.atom_a = atom_a + self.bond_info.atom_b = atom_b + self.bond_info.bond_r = bond_r + + def add_hangle_to_constrain_pair(self, angle_numbers, atom_a, atom_b, atom_c, angle_theta, atom_mass): + """add_hbond_to_constrain_pair""" + self.h_angle_pair = [] + s = 0 + for i in range(angle_numbers): + mass_a = atom_mass[atom_a[i]] + mass_c = atom_mass[atom_c[i]] + if ((0 < mass_a < self.info.constrain_mass) or (0 < mass_c < self.info.constrain_mass)): + for j in range(self.bond_info.bond_numbers): + if ((self.bond_info.atom_a[j] == atom_a[i] and self.bond_info.atom_b[j] == atom_b[i]) + or (self.bond_info.atom_a[j] == atom_b[i] and self.bond_info.atom_b[j] == atom_a[i])): + rab = self.bond_info.bond_r[j] + if ((self.bond_info.atom_a[j] == atom_c[i] and self.bond_info.atom_b[j] == atom_b[i]) + or (self.bond_info.atom_a[j] == atom_b[i] and self.bond_info.atom_b[j] == atom_c[i])): + rbc = self.bond_info.bond_r[j] + + constant_r = math.sqrt(rab * rab + rbc * rbc - 2. * rab * rbc * math.cos(angle_theta[i])) + constrain_k = atom_mass[atom_a[i]] * atom_mass[atom_c[i]] / \ + (atom_mass[atom_a[i]] + atom_mass[atom_c[i]]) + self.h_angle_pair.append([atom_a[i], atom_c[i], constant_r, constrain_k]) + s = s + 1 + self.angle_constrain_pair_numbers = s + + def initial_simple_constrain(self, atom_numbers, dt, box_length, exp_gamma, is_Minimization, atom_mass, + system_freedom): + """initial_simple_constrain""" + self.system_freedom = system_freedom + self.atom_mass = atom_mass + self.info.atom_numbers = atom_numbers + self.info.dt = dt + self.info.quarter_crd_to_uint_crd_cof = [0.25 * self.CONSTANT_UINT_MAX_FLOAT / box_length[0], + 0.25 * self.CONSTANT_UINT_MAX_FLOAT / box_length[1], + 0.25 * self.CONSTANT_UINT_MAX_FLOAT / box_length[2]] + self.info.uint_dr_to_dr_cof = [1 / self.CONSTANT_UINT_MAX_FLOAT * box_length[0], + 1 / self.CONSTANT_UINT_MAX_FLOAT * box_length[1], + 1 / self.CONSTANT_UINT_MAX_FLOAT * box_length[2]] + self.info.volume = box_length[0] * box_length[1] * box_length[2] + self.info.exp_gamma = exp_gamma + self.half_exp_gamma_plus_half = 0.5 * (1. + self.info.exp_gamma) + if is_Minimization: + self.info.exp_gamma = 0.0 + self.info.iteration_numbers = 25 if "iteration_numbers" not in self.controller.Command_Set else int( + self.controller.Command_Set["iteration_numbers"]) + print(" constrain iteration step is ", self.info.iteration_numbers) + self.info.step_length = 1 if "step_length" not in self.controller.Command_Set else float( + self.controller.Command_Set["step_length"]) + print(" constrain step_length is ", self.info.step_length) + self.extra_numbers = 0 + name = self.module_name + "_in_file" + if name in self.controller.Command_Set: + path = self.controller.Command_Set[name] + file = open(path, 'r') + context = file.readlines() + self.extra_numbers = int(context[0].strip()) + file.close() + + self.dt_inverse = 1 / self.info.dt + self.constrain_pair_numbers = self.bond_constrain_pair_numbers + self.angle_constrain_pair_numbers + \ + self.extra_numbers + self.system_freedom -= self.constrain_pair_numbers + + print(" constrain pair number is ", self.constrain_pair_numbers) + + for i in range(self.bond_constrain_pair_numbers): + self.h_constrain_pair.append(self.h_bond_pair[i]) + self.h_constrain_pair[i][-1] = self.info.step_length / self.half_exp_gamma_plus_half \ + * self.h_constrain_pair[i][-1] + for i in range(self.angle_constrain_pair_numbers): + self.h_constrain_pair.append(self.h_bond_pair[i]) + idx = self.bond_constrain_pair_numbers + self.h_constrain_pair[i + idx][-1] = self.info.step_length / \ + self.half_exp_gamma_plus_half * self.h_constrain_pair[i + idx][-1] + if name in self.controller.Command_Set: + path = self.controller.Command_Set[name] + self.read_in_file(path) + + self.is_initialized = 1 + + def read_in_file(self, path): + """read_in_file""" + file = open(path, 'r') + context = file.readlines() + count = self.bond_constrain_pair_numbers + self.angle_constrain_pair_numbers + for i in range(self.extra_numbers): + val = list(map(float, context[i + 1].strip().split())) + atom_i, atom_j, constant_r = int(val[0]), int(val[1]), float(val[2]) + constrain_k = self.info.step_length / self.half_exp_gamma_plus_half * self.atom_mass[atom_i] * \ + self.atom_mass[atom_j] / (self.atom_mass[atom_i] + self.atom_mass[atom_j]) + self.h_constrain_pair.append([atom_i, atom_j, constant_r, constrain_k]) + count += 1 + file.close() diff --git a/model_zoo/research/hpc/sponge/src/simulation.py b/model_zoo/research/hpc/sponge/src/simulation.py index e5474806c61..e02c844c476 100644 --- a/model_zoo/research/hpc/sponge/src/simulation.py +++ b/model_zoo/research/hpc/sponge/src/simulation.py @@ -13,23 +13,29 @@ # limitations under the License. # ============================================================================ '''Simulation''' -import numpy as np -import mindspore.common.dtype as mstype -from mindspore import Tensor -from mindspore import nn -from mindspore.common.parameter import Parameter -from mindspore.ops import functional as F -from mindspore.ops import operations as P +import numpy as np from src.angle import Angle +from src.bd_baro import BD_BARO from src.bond import Bond +from src.crd_molecular_map import CoordinateMolecularMap from src.dihedral import Dihedral from src.langevin_liujian_md import Langevin_Liujian from src.lennard_jones import Lennard_Jones_Information +from src.mc_baro import MC_BARO from src.md_information import md_information from src.nb14 import NON_BOND_14 from src.neighbor_list import neighbor_list from src.particle_mesh_ewald import Particle_Mesh_Ewald +from src.restrain import Restrain_Information +from src.simple_constrain import Simple_Constarin +from src.vatom import Virtual_Information + +import mindspore.common.dtype as mstype +from mindspore import Tensor, nn +from mindspore.common.parameter import Parameter +from mindspore.ops import functional as F +from mindspore.ops import operations as P class controller: @@ -47,6 +53,7 @@ class controller: self.Command_Set = {} self.md_task = None self.commands_from_in_file() + self.punctuation = "," def commands_from_in_file(self): '''command from in file''' @@ -55,10 +62,12 @@ class controller: file.close() self.md_task = context[0].strip() for val in context: - if "=" in val: + val = val.strip() + if val and val[0] != '#' and ("=" in val): + val = val[:val.index(",")] if ',' in val else val assert len(val.strip().split("=")) == 2 flag, value = val.strip().split("=") - value = value.replace(",", '') + value = value.replace(" ", "") flag = flag.replace(" ", "") if flag not in self.Command_Set: self.Command_Set[flag] = value @@ -73,14 +82,99 @@ class Simulation(nn.Cell): super(Simulation, self).__init__() self.control = controller(args_opt) self.md_info = md_information(self.control) - self.bond = Bond(self.control, self.md_info) + self.mode = self.md_info.mode + self.bond = Bond(self.control) + self.bond_is_initialized = self.bond.is_initialized self.angle = Angle(self.control) + self.angle_is_initialized = self.angle.is_initialized self.dihedral = Dihedral(self.control) + self.dihedral_is_initialized = self.dihedral.is_initialized self.nb14 = NON_BOND_14(self.control, self.dihedral, self.md_info.atom_numbers) + self.nb14_is_initialized = self.nb14.is_initialized self.nb_info = neighbor_list(self.control, self.md_info.atom_numbers, self.md_info.box_length) - self.LJ_info = Lennard_Jones_Information(self.control) + self.LJ_info = Lennard_Jones_Information(self.control, self.md_info.nb.cutoff, self.md_info.sys.box_length) + self.LJ_info_is_initialized = self.LJ_info.is_initialized + self.liujian_info = Langevin_Liujian(self.control, self.md_info.atom_numbers) + self.liujian_info_is_initialized = self.liujian_info.is_initialized self.pme_method = Particle_Mesh_Ewald(self.control, self.md_info) + self.pme_is_initialized = self.pme_method.is_initialized + self.restrain = Restrain_Information(self.control, self.md_info.atom_numbers, self.md_info.crd) + self.restrain_is_initialized = self.restrain.is_initialized + self.simple_constrain_is_initialized = 0 + + self.simple_constrain = Simple_Constarin(self.control, self.md_info, self.bond, self.angle, self.liujian_info) + self.simple_constrain_is_initialized = self.simple_constrain.is_initialized + self.freedom = self.simple_constrain.system_freedom + + self.vatom = Virtual_Information(self.control, self.md_info, self.md_info.sys.freedom) + self.vatom_is_initialized = 1 + + self.random = P.UniformReal(seed=1) + self.pow = P.Pow() + + self.mol_map = CoordinateMolecularMap(self.md_info.atom_numbers, self.md_info.sys.box_length, self.md_info.crd, + self.md_info.nb.excluded_atom_numbers, self.md_info.nb.h_excluded_numbers, + self.md_info.nb.h_excluded_list_start, self.md_info.nb.h_excluded_list) + self.mol_map_is_initialized = 1 + self.init_params() + self.init_Tensor() + self.op_define() + self.op_define_2() + self.depend = P.Depend() + self.print = P.Print() + self.total_count = Parameter(Tensor(0, mstype.int32), requires_grad=False) + self.accept_count = Parameter(Tensor(0, mstype.int32), requires_grad=False) + self.is_molecule_map_output = self.md_info.output.is_molecule_map_output + self.target_pressure = self.md_info.sys.target_pressure + self.Nx = self.nb_info.Nx + self.Ny = self.nb_info.Ny + self.Nz = self.nb_info.Nz + self.PME_inverse_box_vector = Parameter(Tensor(self.pme_method.PME_inverse_box_vector, mstype.float32), + requires_grad=False) + self.mc_baro_is_initialized = 0 + self.bd_baro_is_initialized = 0 + + if self.mode == 2 and self.control.Command_Set["barostat"] == "monte_carlo": + self.mc_baro = MC_BARO(self.control, self.md_info.atom_numbers, self.md_info.sys.target_pressure, + self.md_info.sys.box_length, self.md_info.res.is_initialized, self.md_info.mode) + self.mc_baro_is_initialized = self.mc_baro.is_initialized + self.update_interval = self.mc_baro.update_interval + self.mc_baro_energy_old = Parameter(Tensor(0, mstype.float32), requires_grad=False) + self.potential = Parameter(Tensor(0, mstype.float32), requires_grad=False) + self.frc_backup = Parameter(Tensor(np.zeros([self.atom_numbers, 3]), mstype.float32), requires_grad=False) + self.crd_backup = Parameter(Tensor(np.zeros([self.atom_numbers, 3]), mstype.float32), requires_grad=False) + self.crd_scale_factor = Parameter(Tensor(0.0, mstype.float32), requires_grad=False) + self.system_reinitializing_count = Parameter(Tensor(0, mstype.int32), requires_grad=False) + self.mc_baro_energy_new = Parameter(Tensor(0.0, mstype.float32), requires_grad=False) + self.scale_coordinate_by_residue = Parameter(Tensor(0, mstype.float32), requires_grad=False) + self.extra_term = Parameter(Tensor(0, mstype.float32), requires_grad=False) + self.DeltaV = Parameter(Tensor(0.0, mstype.float32), requires_grad=False) + self.target_temperature = self.md_info.sys.target_temperature + self.VDevided = Parameter(Tensor(0.0, mstype.float32), requires_grad=False) + self.log = P.Log() + self.mc_baro_accept_possibility = Parameter(Tensor(0, mstype.float32), requires_grad=False) + self.exp = P.Exp() + self.mc_baro_newV = self.mc_baro.newV + self.mc_baro_V0 = Parameter(Tensor(self.mc_baro.V0, mstype.float32), requires_grad=False) + self.mc_baro_newV = self.mc_baro.newV + self.check_interval = self.mc_baro.check_interval + + if self.mode == 2 and self.control.Command_Set["barostat"] == "berendsen": + self.bd_baro = BD_BARO(self.control, self.md_info.sys.target_pressure, self.md_info.sys.box_length, + self.md_info.mode) + self.bd_baro_is_initialized = self.bd_baro.is_initialized + self.update_interval = self.bd_baro.update_interval + self.pressure = Parameter(Tensor(self.md_info.sys.d_pressure, mstype.float32), requires_grad=False) + self.compressibility = self.bd_baro.compressibility + self.bd_baro_dt = self.bd_baro.dt + self.bd_baro_taup = self.bd_baro.taup + self.system_reinitializing_count = Parameter(Tensor(0, mstype.int32), requires_grad=False) + self.bd_baro_newV = Parameter(Tensor(self.bd_baro.newV, mstype.float32), requires_grad=False) + self.bd_baro_V0 = Parameter(Tensor(self.bd_baro.V0, mstype.float32), requires_grad=False) + + def init_params(self): + """init_params""" self.bond_energy_sum = Tensor(0, mstype.int32) self.angle_energy_sum = Tensor(0, mstype.int32) self.dihedral_energy_sum = Tensor(0, mstype.int32) @@ -101,7 +195,8 @@ class Simulation(nn.Cell): self.grid_numbers = self.nb_info.grid_numbers self.max_atom_in_grid_numbers = self.nb_info.max_atom_in_grid_numbers self.max_neighbor_numbers = self.nb_info.max_neighbor_numbers - self.excluded_atom_numbers = self.nb_info.excluded_atom_numbers + # self.excluded_atom_numbers = self.nb_info.excluded_atom_numbers + self.excluded_atom_numbers = self.md_info.nb.excluded_atom_numbers self.refresh_count = Parameter(Tensor(self.nb_info.refresh_count, mstype.int32), requires_grad=False) self.refresh_interval = self.nb_info.refresh_interval self.skin = self.nb_info.skin @@ -115,24 +210,39 @@ class Simulation(nn.Cell): self.fftx = self.pme_method.fftx self.ffty = self.pme_method.ffty self.fftz = self.pme_method.fftz - self.random_seed = self.liujian_info.rand_seed + self.random_seed = self.liujian_info.random_seed self.dt = self.liujian_info.dt self.half_dt = self.liujian_info.half_dt self.exp_gamma = self.liujian_info.exp_gamma - self.init_Tensor() - self.op_define() self.update = False self.file = None self.datfile = None + self.max_velocity = self.liujian_info.max_velocity + + # bingshui + self.CONSTANT_kB = 0.00198716 def init_Tensor(self): '''init tensor''' + # MD_Reset_Atom_Energy_And_Virial + self.uint_crd = Parameter(Tensor(np.zeros([self.atom_numbers, 3], dtype=np.uint32), mstype.uint32), + requires_grad=False) + self.need_potential = Tensor(0, mstype.int32) + self.need_pressure = Tensor(0, mstype.int32) + # self.potential = Tensor(0, mstype.float32) + self.atom_energy = Parameter(Tensor([0] * self.atom_numbers, mstype.float32), requires_grad=False) + self.atom_virial = Parameter(Tensor([0] * self.atom_numbers, mstype.float32), requires_grad=False) + self.frc = Parameter(Tensor(np.zeros([self.atom_numbers, 3]), mstype.float32), requires_grad=False) + self.crd = Parameter( - Tensor(np.float32(np.asarray(self.md_info.coordinate).reshape([self.atom_numbers, 3])), mstype.float32), + Tensor(np.array(self.md_info.coordinate).reshape([self.atom_numbers, 3]), mstype.float32), requires_grad=False) - self.crd_to_uint_crd_cof = Tensor(np.asarray(self.md_info.crd_to_uint_crd_cof, np.float32), mstype.float32) - self.uint_dr_to_dr_cof = Parameter( - Tensor(np.asarray(self.md_info.uint_dr_to_dr_cof, np.float32), mstype.float32), requires_grad=False) + self.crd_to_uint_crd_cof = Tensor(np.asarray(self.md_info.pbc.crd_to_uint_crd_cof, np.float32), mstype.float32) + self.quarter_crd_to_uint_crd_cof = Tensor(np.asarray(self.md_info.pbc.quarter_crd_to_uint_crd_cof, np.float32), + mstype.float32) + + self.uint_dr_to_dr_cof = Parameter(Tensor(self.md_info.pbc.uint_dr_to_dr_cof, mstype.float32), + requires_grad=False) self.box_length = Tensor(self.md_info.box_length, mstype.float32) self.charge = Parameter(Tensor(np.asarray(self.md_info.h_charge, dtype=np.float32), mstype.float32), requires_grad=False) @@ -140,12 +250,13 @@ class Simulation(nn.Cell): requires_grad=False) self.last_crd = Parameter(Tensor(np.zeros([self.atom_numbers, 3], dtype=np.float32), mstype.float32), requires_grad=False) - self.uint_crd = Parameter(Tensor(np.zeros([self.atom_numbers, 3], dtype=np.uint32), mstype.uint32), - requires_grad=False) + self.mass = Tensor(self.md_info.h_mass, mstype.float32) self.mass_inverse = Tensor(self.md_info.h_mass_inverse, mstype.float32) + self.res_mass = Tensor(self.md_info.res.h_mass, mstype.float32) + self.res_mass_inverse = Tensor(self.md_info.res.h_mass_inverse, mstype.float32) + self.res_start = Tensor(self.md_info.h_res_start, mstype.int32) self.res_end = Tensor(self.md_info.h_res_end, mstype.int32) - self.mass = Tensor(self.md_info.h_mass, mstype.float32) self.velocity = Parameter(Tensor(self.md_info.velocity, mstype.float32), requires_grad=False) self.acc = Parameter(Tensor(np.zeros([self.atom_numbers, 3], np.float32), mstype.float32), requires_grad=False) self.bond_atom_a = Tensor(np.asarray(self.bond.h_atom_a, np.int32), mstype.int32) @@ -161,17 +272,19 @@ class Simulation(nn.Cell): self.dihedral_atom_b = Tensor(np.asarray(self.dihedral.h_atom_b, np.int32), mstype.int32) self.dihedral_atom_c = Tensor(np.asarray(self.dihedral.h_atom_c, np.int32), mstype.int32) self.dihedral_atom_d = Tensor(np.asarray(self.dihedral.h_atom_d, np.int32), mstype.int32) - self.pk = Tensor(np.asarray(self.dihedral.pk, np.float32), mstype.float32) - self.gamc = Tensor(np.asarray(self.dihedral.gamc, np.float32), mstype.float32) - self.gams = Tensor(np.asarray(self.dihedral.gams, np.float32), mstype.float32) - self.pn = Tensor(np.asarray(self.dihedral.pn, np.float32), mstype.float32) - self.ipn = Tensor(np.asarray(self.dihedral.ipn, np.int32), mstype.int32) + self.pk = Tensor(np.asarray(self.dihedral.h_pk, np.float32), mstype.float32) + self.gamc = Tensor(np.asarray(self.dihedral.h_gamc, np.float32), mstype.float32) + self.gams = Tensor(np.asarray(self.dihedral.h_gams, np.float32), mstype.float32) + self.pn = Tensor(np.asarray(self.dihedral.h_pn, np.float32), mstype.float32) + self.ipn = Tensor(np.asarray(self.dihedral.h_ipn, np.int32), mstype.int32) self.nb14_atom_a = Tensor(np.asarray(self.nb14.h_atom_a, np.int32), mstype.int32) self.nb14_atom_b = Tensor(np.asarray(self.nb14.h_atom_b, np.int32), mstype.int32) self.lj_scale_factor = Tensor(np.asarray(self.nb14.h_lj_scale_factor, np.float32), mstype.float32) self.cf_scale_factor = Tensor(np.asarray(self.nb14.h_cf_scale_factor, np.float32), mstype.float32) self.grid_N = Tensor(self.nb_info.grid_N, mstype.int32) - self.grid_length_inverse = Tensor(self.nb_info.grid_length_inverse, mstype.float32) + self.grid_length = Parameter(Tensor(self.nb_info.grid_length, mstype.float32), requires_grad=False) + self.grid_length_inverse = Parameter(Tensor(self.nb_info.grid_length_inverse, mstype.float32), + requires_grad=False) self.bucket = Parameter(Tensor( np.asarray(self.nb_info.bucket, np.int32).reshape([self.grid_numbers, self.max_atom_in_grid_numbers]), mstype.int32), requires_grad=False) @@ -187,24 +300,29 @@ class Simulation(nn.Cell): self.nl_atom_serial = Parameter( Tensor(np.zeros([self.atom_numbers, self.max_neighbor_numbers], np.int32), mstype.int32), requires_grad=False) - self.excluded_list_start = Tensor(np.asarray(self.nb_info.excluded_list_start, np.int32), mstype.int32) - self.excluded_list = Tensor(np.asarray(self.nb_info.excluded_list, np.int32), mstype.int32) - self.excluded_numbers = Tensor(np.asarray(self.nb_info.excluded_numbers, np.int32), mstype.int32) + self.excluded_list_start = Tensor(np.asarray(self.md_info.nb.h_excluded_list_start, np.int32), mstype.int32) + self.excluded_list = Tensor(np.asarray(self.md_info.nb.h_excluded_list, np.int32), mstype.int32) + self.excluded_numbers = Tensor(np.asarray(self.md_info.nb.h_excluded_numbers, np.int32), mstype.int32) + self.need_refresh_flag = Tensor(np.asarray([0], np.int32), mstype.int32) - self.atom_LJ_type = Tensor(np.asarray(self.LJ_info.atom_LJ_type, dtype=np.int32), mstype.int32) - self.LJ_A = Tensor(np.asarray(self.LJ_info.LJ_A, dtype=np.float32), mstype.float32) - self.LJ_B = Tensor(np.asarray(self.LJ_info.LJ_B, dtype=np.float32), mstype.float32) + self.atom_LJ_type = Tensor(self.LJ_info.atom_LJ_type, mstype.int32) + self.LJ_A = Tensor(self.LJ_info.h_LJ_A, mstype.float32) + self.LJ_B = Tensor(self.LJ_info.h_LJ_B, mstype.float32) self.sqrt_mass = Tensor(self.liujian_info.h_sqrt_mass, mstype.float32) self.rand_state = Parameter(Tensor(self.liujian_info.rand_state, mstype.float32)) self.zero_fp_tensor = Tensor(np.asarray([0,], np.float32)) + self.zero_frc = Parameter(Tensor(np.zeros([self.atom_numbers, 3], dtype=np.float32), mstype.float32), + requires_grad=False) def op_define(self): '''op define''' self.crd_to_uint_crd = P.CrdToUintCrd(self.atom_numbers) + self.crd_to_uint_crd_quarter = P.CrdToUintCrdQuarter(self.atom_numbers) self.mdtemp = P.MDTemperature(self.residue_numbers, self.atom_numbers) self.setup_random_state = P.MDIterationSetupRandState(self.atom_numbers, self.random_seed) - self.bond_force_with_atom_energy = P.BondForceWithAtomEnergy(bond_numbers=self.bond_numbers, - atom_numbers=self.atom_numbers) + + self.bond_force_with_atom_energy_virial = P.BondForceWithAtomEnergyAndVirial(bond_numbers=self.bond_numbers, + atom_numbers=self.atom_numbers) self.angle_force_with_atom_energy = P.AngleForceWithAtomEnergy(angle_numbers=self.angle_numbers) self.dihedral_force_with_atom_energy = P.DihedralForceWithAtomEnergy(dihedral_numbers=self.dihedral_numbers) self.nb14_force_with_atom_energy = P.Dihedral14LJCFForceWithAtomEnergy(nb14_numbers=self.nb14_numbers, @@ -215,7 +333,6 @@ class Simulation(nn.Cell): self.pme_reciprocal_force = P.PMEReciprocalForce(self.atom_numbers, self.beta, self.fftx, self.ffty, self.fftz, self.md_info.box_length[0], self.md_info.box_length[1], self.md_info.box_length[2]) - self.bond_energy = P.BondEnergy(self.bond_numbers, self.atom_numbers) self.angle_energy = P.AngleEnergy(self.angle_numbers) self.dihedral_energy = P.DihedralEnergy(self.dihedral_numbers) @@ -225,77 +342,204 @@ class Simulation(nn.Cell): self.pme_energy = P.PMEEnergy(self.atom_numbers, self.excluded_atom_numbers, self.beta, self.fftx, self.ffty, self.fftz, self.md_info.box_length[0], self.md_info.box_length[1], self.md_info.box_length[2]) - self.md_iteration_leap_frog_liujian = P.MDIterationLeapFrogLiujian(self.atom_numbers, self.half_dt, self.dt, self.exp_gamma) - self.neighbor_list_update_init = P.NeighborListUpdate(grid_numbers=self.grid_numbers, - atom_numbers=self.atom_numbers, not_first_time=0, - nxy=self.nxy, - excluded_atom_numbers=self.excluded_atom_numbers, - cutoff_square=self.cutoff_square, - half_skin_square=self.half_skin_square, - cutoff_with_skin=self.cutoff_with_skin, - half_cutoff_with_skin=self.half_cutoff_with_skin, - cutoff_with_skin_square=self.cutoff_with_skin_square, - refresh_interval=self.refresh_interval, - cutoff=self.cutoff, skin=self.skin, - max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, - max_neighbor_numbers=self.max_neighbor_numbers) + self.md_iteration_leap_frog_liujian_with_max_vel = P.MDIterationLeapFrogLiujianWithMaxVel(self.atom_numbers, + self.half_dt, self.dt, + self.exp_gamma, + self.max_velocity) + self.neighbor_list_update = \ + P.NeighborListUpdate(grid_numbers=self.grid_numbers, + atom_numbers=self.atom_numbers, + not_first_time=1, nxy=self.nxy, + excluded_atom_numbers=self.excluded_atom_numbers, + cutoff_square=self.cutoff_square, + half_skin_square=self.half_skin_square, + cutoff_with_skin=self.cutoff_with_skin, + half_cutoff_with_skin=self.half_cutoff_with_skin, + cutoff_with_skin_square=self.cutoff_with_skin_square, + refresh_interval=self.refresh_interval, cutoff=self.cutoff, + skin=self.skin, + max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, + max_neighbor_numbers=self.max_neighbor_numbers) + + self.neighbor_list_update_forced_update = \ + P.NeighborListUpdate(grid_numbers=self.grid_numbers, + atom_numbers=self.atom_numbers, + not_first_time=1, nxy=self.nxy, + excluded_atom_numbers=self.excluded_atom_numbers, + cutoff_square=self.cutoff_square, + half_skin_square=self.half_skin_square, + cutoff_with_skin=self.cutoff_with_skin, + half_cutoff_with_skin=self.half_cutoff_with_skin, + cutoff_with_skin_square=self.cutoff_with_skin_square, + refresh_interval=self.refresh_interval, + cutoff=self.cutoff, + skin=self.skin, + max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, + max_neighbor_numbers=self.max_neighbor_numbers, + forced_update=1) + + self.neighbor_list_update_nb = \ + P.NeighborListUpdate(grid_numbers=self.grid_numbers, + atom_numbers=self.atom_numbers, + not_first_time=1, nxy=self.nxy, + excluded_atom_numbers=self.excluded_atom_numbers, + cutoff_square=self.cutoff_square, + half_skin_square=self.half_skin_square, + cutoff_with_skin=self.cutoff_with_skin, + half_cutoff_with_skin=self.half_cutoff_with_skin, + cutoff_with_skin_square=self.cutoff_with_skin_square, + refresh_interval=self.refresh_interval, + cutoff=self.cutoff, + skin=self.skin, + max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, + max_neighbor_numbers=self.max_neighbor_numbers, + forced_update=1, forced_check=1) + + def op_define_2(self): + """op_define_2""" + self.neighbor_list_update_mc = P.NeighborListUpdate(grid_numbers=self.grid_numbers, + atom_numbers=self.atom_numbers, + not_first_time=1, nxy=self.nxy, + excluded_atom_numbers=self.excluded_atom_numbers, + cutoff_square=self.cutoff_square, + half_skin_square=self.half_skin_square, + cutoff_with_skin=self.cutoff_with_skin, + half_cutoff_with_skin=self.half_cutoff_with_skin, + cutoff_with_skin_square=self.cutoff_with_skin_square, + refresh_interval=self.refresh_interval, + cutoff=self.cutoff, + skin=self.skin, + max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, + max_neighbor_numbers=self.max_neighbor_numbers, + forced_update=0, forced_check=1) - self.neighbor_list_update = P.NeighborListUpdate(grid_numbers=self.grid_numbers, atom_numbers=self.atom_numbers, - not_first_time=1, nxy=self.nxy, - excluded_atom_numbers=self.excluded_atom_numbers, - cutoff_square=self.cutoff_square, - half_skin_square=self.half_skin_square, - cutoff_with_skin=self.cutoff_with_skin, - half_cutoff_with_skin=self.half_cutoff_with_skin, - cutoff_with_skin_square=self.cutoff_with_skin_square, - refresh_interval=self.refresh_interval, cutoff=self.cutoff, - skin=self.skin, - max_atom_in_grid_numbers=self.max_atom_in_grid_numbers, - max_neighbor_numbers=self.max_neighbor_numbers) self.random_force = Tensor(np.zeros([self.atom_numbers, 3], np.float32), mstype.float32) + # simple_constrain + self.constrain_pair_numbers = self.simple_constrain.constrain_pair_numbers + self.last_pair_dr = Parameter(Tensor(np.zeros([self.constrain_pair_numbers, 3], np.float32), mstype.float32), + requires_grad=False) + if self.simple_constrain_is_initialized: + self.constrain_pair_numbers = self.simple_constrain.constrain_pair_numbers + self.last_crd_to_dr = P.lastcrdtodr(self.atom_numbers, self.constrain_pair_numbers) + self.constrain_pair = np.array(self.simple_constrain.h_constrain_pair) + self.atom_i_serials = Tensor(self.constrain_pair[:, 0], mstype.int32) + self.atom_j_serials = Tensor(self.constrain_pair[:, 1], mstype.int32) + self.constant_rs = Tensor(self.constrain_pair[:, 2], mstype.float32) + self.constrain_ks = Tensor(self.constrain_pair[:, 3], mstype.float32) + self.last_pair_dr = Parameter( + Tensor(np.zeros([self.constrain_pair_numbers, 3], np.float32), mstype.float32), requires_grad=False) + self.constrain_frc = Parameter(Tensor(np.zeros([self.atom_numbers, 3], np.float32), mstype.float32), + requires_grad=False) + self.iteration_numbers = self.simple_constrain.info.iteration_numbers + self.half_exp_gamma_plus_half = self.simple_constrain.half_exp_gamma_plus_half + self.refresh_uint_crd = P.refreshuintcrd(self.atom_numbers, self.half_exp_gamma_plus_half) + self.need_pressure = 0 + self.constrain_force_cycle_with_virial = P.constrainforcecyclewithvirial(self.atom_numbers, + self.constrain_pair_numbers) + self.constrain_force_cycle = P.ConstrainForceCycle(self.atom_numbers, self.constrain_pair_numbers) + self.dt_inverse = self.simple_constrain.dt_inverse + self.refresh_crd_vel = P.refreshcrdvel(self.atom_numbers, self.dt_inverse, self.dt, self.exp_gamma, + self.half_exp_gamma_plus_half) + + if self.mol_map_is_initialized: + self.refresh_boxmaptimes = P.refreshboxmaptimes(self.atom_numbers) + self.box_map_times = Parameter(Tensor(self.mol_map.h_box_map_times, mstype.int32), requires_grad=False) + self.residue_numbers = self.md_info.residue_numbers + self.getcenterofmass = P.GetCenterOfMass(self.residue_numbers) + self.mapcenterofmass = P.MapCenterOfMass(self.residue_numbers, scaler=1.0) + + self.md_iteration_leap_frog = P.MDIterationLeapFrog(self.atom_numbers, self.dt) + self.md_iteration_leap_frog_with_max_vel = P.MDIterationLeapFrogWithMaxVel(self.atom_numbers, self.dt, + self.max_velocity) + self.md_information_gradient_descent = P.MDIterationGradientDescent(self.atom_numbers, self.dt * self.dt) + def Simulation_Beforce_Caculate_Force(self): '''simulation before calculate force''' - crd_to_uint_crd_cof = 0.5 * self.crd_to_uint_crd_cof - uint_crd = self.crd_to_uint_crd(crd_to_uint_crd_cof, self.crd) - return uint_crd + self.uint_crd = self.crd_to_uint_crd_quarter(self.quarter_crd_to_uint_crd_cof, self.crd) + return self.uint_crd def Simulation_Caculate_Force(self, uint_crd, scaler, nl_atom_numbers, nl_atom_serial): '''simulation calculate force''' - bond_force, _ = self.bond_force_with_atom_energy(uint_crd, scaler, self.bond_atom_a, - self.bond_atom_b, self.bond_k, self.bond_r0) + uint_crd = self.Simulation_Beforce_Caculate_Force() + force = self.zero_frc + if self.LJ_info_is_initialized: + lj_force = self.lj_force_pme_direct_force(uint_crd, self.atom_LJ_type, self.charge, scaler, nl_atom_numbers, + nl_atom_serial, self.LJ_A, self.LJ_B) + force = force + lj_force - angle_force, _ = self.angle_force_with_atom_energy(uint_crd, scaler, self.angle_atom_a, - self.angle_atom_b, self.angle_atom_c, - self.angle_k, self.angle_theta0) + if self.pme_is_initialized: + pme_excluded_force = self.pme_excluded_force(uint_crd, scaler, self.charge, self.excluded_list_start, + self.excluded_list, self.excluded_numbers) - dihedral_force, _ = self.dihedral_force_with_atom_energy(uint_crd, scaler, - self.dihedral_atom_a, - self.dihedral_atom_b, - self.dihedral_atom_c, - self.dihedral_atom_d, self.ipn, - self.pk, self.gamc, self.gams, - self.pn) + pme_reciprocal_force = self.pme_reciprocal_force(uint_crd, self.charge) + force = force + pme_excluded_force + pme_reciprocal_force + if self.nb14_is_initialized: + nb14_force, _ = self.nb14_force_with_atom_energy(uint_crd, self.atom_LJ_type, self.charge, + scaler, self.nb14_atom_a, self.nb14_atom_b, + self.lj_scale_factor, self.cf_scale_factor, + self.LJ_A, self.LJ_B) + force = force + nb14_force - nb14_force, _ = self.nb14_force_with_atom_energy(uint_crd, self.atom_LJ_type, self.charge, - scaler, self.nb14_atom_a, self.nb14_atom_b, - self.lj_scale_factor, self.cf_scale_factor, - self.LJ_A, self.LJ_B) + if self.bond_is_initialized: + bond_force, _, _ = self.bond_force_with_atom_energy_virial(uint_crd, scaler, self.bond_atom_a, + self.bond_atom_b, self.bond_k, self.bond_r0) + force = force + bond_force + if self.angle_is_initialized: + angle_force, _ = self.angle_force_with_atom_energy(uint_crd, scaler, self.angle_atom_a, + self.angle_atom_b, self.angle_atom_c, + self.angle_k, self.angle_theta0) + force = force + angle_force + if self.dihedral_is_initialized: + dihedral_force, _ = self.dihedral_force_with_atom_energy(uint_crd, scaler, + self.dihedral_atom_a, + self.dihedral_atom_b, + self.dihedral_atom_c, + self.dihedral_atom_d, self.ipn, + self.pk, self.gamc, self.gams, + self.pn) + force = force + dihedral_force + + if self.restrain_is_initialized: + _, _, restrain_frc = self.restrain_force_with_atom_energy_and_virial(self.restrain_list, + self.crd, + self.crd_ref, + self.box_length) + force = force + restrain_frc - lj_force = self.lj_force_pme_direct_force(uint_crd, self.atom_LJ_type, self.charge, scaler, nl_atom_numbers, - nl_atom_serial, self.LJ_A, self.LJ_B) - pme_excluded_force = self.pme_excluded_force(uint_crd, scaler, self.charge, self.excluded_list_start, - self.excluded_list, self.excluded_numbers) - pme_reciprocal_force = self.pme_reciprocal_force(uint_crd, self.charge) - force = P.AddN()( - [bond_force, angle_force, dihedral_force, nb14_force, lj_force, pme_excluded_force, pme_reciprocal_force]) return force def Simulation_Caculate_Energy(self, uint_crd, uint_dr_to_dr_cof): '''simulation calculate energy''' + + lj_energy = self.lj_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, self.nl_atom_numbers, + self.nl_atom_serial, self.LJ_A, self.LJ_B) + + lj_energy_sum = P.ReduceSum(True)(lj_energy) + # lj_energy_sum = self.zero_fp_tensor + + reciprocal_energy, self_energy, direct_energy, correction_energy = self.pme_energy(uint_crd, self.charge, + self.nl_atom_numbers, + self.nl_atom_serial, + uint_dr_to_dr_cof, + self.excluded_list_start, + self.excluded_list, + self.excluded_numbers) + ee_ene = reciprocal_energy + self_energy + direct_energy + correction_energy + # ee_ene = self.zero_fp_tensor + + nb14_lj_energy = self.nb14_lj_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, + self.nb14_atom_a, self.nb14_atom_b, self.lj_scale_factor, self.LJ_A, + self.LJ_B) + nb14_cf_energy = self.nb14_cf_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, + self.nb14_atom_a, self.nb14_atom_b, self.cf_scale_factor) + nb14_lj_energy_sum = P.ReduceSum(True)(nb14_lj_energy) + nb14_cf_energy_sum = P.ReduceSum(True)(nb14_cf_energy) + # nb14_lj_energy_sum = self.zero_fp_tensor + # nb14_cf_energy_sum = self.zero_fp_tensor bond_energy = self.bond_energy(uint_crd, uint_dr_to_dr_cof, self.bond_atom_a, self.bond_atom_b, self.bond_k, self.bond_r0) bond_energy_sum = P.ReduceSum(True)(bond_energy) @@ -309,26 +553,6 @@ class Simulation(nn.Cell): self.gams, self.pn) dihedral_energy_sum = P.ReduceSum(True)(dihedral_energy) - nb14_lj_energy = self.nb14_lj_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, - self.nb14_atom_a, self.nb14_atom_b, self.lj_scale_factor, self.LJ_A, - self.LJ_B) - nb14_cf_energy = self.nb14_cf_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, - self.nb14_atom_a, self.nb14_atom_b, self.cf_scale_factor) - nb14_lj_energy_sum = P.ReduceSum(True)(nb14_lj_energy) - nb14_cf_energy_sum = P.ReduceSum(True)(nb14_cf_energy) - - lj_energy = self.lj_energy(uint_crd, self.atom_LJ_type, self.charge, uint_dr_to_dr_cof, self.nl_atom_numbers, - self.nl_atom_serial, self.LJ_A, self.LJ_B) - lj_energy_sum = P.ReduceSum(True)(lj_energy) - - reciprocal_energy, self_energy, direct_energy, correction_energy = self.pme_energy(uint_crd, self.charge, - self.nl_atom_numbers, - self.nl_atom_serial, - uint_dr_to_dr_cof, - self.excluded_list_start, - self.excluded_list, - self.excluded_numbers) - ee_ene = reciprocal_energy + self_energy + direct_energy + correction_energy total_energy = P.AddN()( [bond_energy_sum, angle_energy_sum, dihedral_energy_sum, nb14_lj_energy_sum, nb14_cf_energy_sum, lj_energy_sum, ee_ene]) @@ -336,19 +560,43 @@ class Simulation(nn.Cell): lj_energy_sum, ee_ene, total_energy def Simulation_Temperature(self): - '''caculate temperature''' + """calculate temperature""" res_ek_energy = self.mdtemp(self.res_start, self.res_end, self.velocity, self.mass) temperature = P.ReduceSum()(res_ek_energy) return temperature def Simulation_MDIterationLeapFrog_Liujian(self, inverse_mass, sqrt_mass_inverse, crd, frc, rand_state, random_frc): '''simulation leap frog iteration liujian''' - crd = self.md_iteration_leap_frog_liujian(inverse_mass, sqrt_mass_inverse, self.velocity, crd, frc, self.acc, - rand_state, random_frc) + if self.max_velocity <= 0: + crd = self.md_iteration_leap_frog_liujian(inverse_mass, sqrt_mass_inverse, self.velocity, crd, frc, + self.acc, + rand_state, random_frc) + else: + crd = self.md_iteration_leap_frog_liujian_with_max_vel(inverse_mass, sqrt_mass_inverse, self.velocity, crd, + frc, self.acc, + rand_state, random_frc) vel = F.depend(self.velocity, crd) acc = F.depend(self.acc, crd) return vel, crd, acc + def Simulation_MDIterationLeapFrog(self, force): + '''simulation leap frog''' + if self.max_velocity <= 0: + res = self.md_iteration_leap_frog(self.velocity, self.crd, force, self.acc, self.mass_inverse) + else: + res = self.md_iteration_leap_frog_with_max_vel(self.velocity, self.crd, force, self.acc, self.mass_inverse) + vel = F.depend(self.velocity, res) + crd = F.depend(self.crd, res) + return vel, crd, res + + def Simulation_MDInformationGradientDescent(self, force): + # print("Simulation_MDInformationGradientDescent") + res = self.md_information_gradient_descent(self.crd, force) + self.velocity = self.zero_frc + vel = F.depend(self.velocity, res) + crd = F.depend(self.crd, res) + return vel, crd, res + def Main_Print(self, *args): """compute the temperature""" steps, temperature, total_potential_energy, sigma_of_bond_ene, sigma_of_angle_ene, sigma_of_dihedral_ene, \ @@ -359,7 +607,7 @@ class Simulation(nn.Cell): temperature = temperature.asnumpy() total_potential_energy = total_potential_energy.asnumpy() - print("{:>7.0f} {:>7.3f} {:>11.3f}".format(steps, float(temperature), float(total_potential_energy)), + print("{:>7.0f} {:>7.3f} {:>11.3f}".format(steps + 1, float(temperature), float(total_potential_energy)), end=" ") if self.bond.bond_numbers > 0: sigma_of_bond_ene = sigma_of_bond_ene.asnumpy() @@ -405,34 +653,304 @@ class Simulation(nn.Cell): self.datfile.close() print("Save .dat file successfully!") + # 控压部分代码 + def Volume_Change_Attempt(self, boxlength, DeltaV_max): + """Volume_Change_Attempt""" + nrand = self.random((1, 1)) + DeltaV = nrand * DeltaV_max + V = boxlength[0] * boxlength[1] * boxlength[2] + # crd_scale_factor = Tensor(np.crbt((V + DeltaV) / V), mstype.float32) + crd_scale_factor = self.pow((V + DeltaV) / V, -3) + return crd_scale_factor + + def Update_Volume(self, factor): + """Update_Volume""" + self.CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + # f_inv = 1.0 / factor + self.box_length = factor * self.box_length + self.crd_to_uint_crd_cof = self.CONSTANT_UINT_MAX_FLOAT / self.box_length + self.quarter_crd_to_uint_crd_cof = 0.25 * self.crd_to_uint_crd_cof + self.uint_dr_to_dr_cof = 1.0 / self.crd_to_uint_crd_cof + self.uint_crd = self.crd_to_uint_crd_quarter(self.quarter_crd_to_uint_crd_cof, self.crd) + + def Neighbor_List_Update_Volume(self, box_length): + """Neighbor_List_Update_Volume""" + self.quarter_crd_to_uint_crd_cof = 0.25 * self.CONSTANT_UINT_MAX_FLOAT / box_length + self.uint_dr_to_dr_cof = 1.0 / self.CONSTANT_UINT_MAX_FLOAT * box_length + self.grid_length[0] = box_length[0] / self.Nx + self.grid_length[1] = box_length[1] / self.Ny + self.grid_length[2] = box_length[1] / self.Nz + self.grid_length_inverse = 1.0 / self.grid_length + + def LJ_Update_Volume(self): + """main destroy""" + if self.LJ_info_is_initialized: + # self.uint_dr_to_dr_cof = 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length + self.volume = self.box_length[0] * self.box_length[1] * self.box_length[2] + + def PME_Update_Volume(self, factor): + """PME_Update_Volume""" + factor_inverse = 1.0 / factor + self.PME_inverse_box_vector[0] = self.fftx / self.box_length[0] + self.PME_inverse_box_vector[1] = self.ffty / self.box_length[1] + self.PME_inverse_box_vector[2] = self.fftz / self.box_length[2] + self.PME_inverse_box_vector = factor_inverse * self.PME_inverse_box_vector + self.beta = self.beta * factor + # self.PME_BC = self.PME_BC * factor_inverse #scale list + self.neutralizing_factor = self.pow(factor, 5.0) + + def Simple_Constrain_Update_Volume(self): + """Simple_Constrain_Update_Volume""" + if self.simple_constrain_is_initialized: + self.quarter_crd_to_uint_crd_cof = 0.25 * self.CONSTANT_UINT_MAX_FLOAT / self.box_length + self.uint_dr_to_dr_cof = 1.0 / self.CONSTANT_UINT_MAX_FLOAT * self.box_length + self.volume = self.box_length[0] * self.box_length[1] * self.box_length[2] + + def Main_Volume_Change(self, factor): + """Main_Volume_Change""" + self.Update_Volume(factor) + self.Neighbor_List_Update_Volume(self.box_length) + _ = self.neighbor_list_update_nb(self.atom_numbers_in_grid_bucket, self.bucket, + self.crd, self.box_length, self.grid_N, + self.grid_length_inverse, self.atom_in_grid_serial, + self.old_crd, self.crd_to_uint_crd_cof, self.uint_crd, + self.pointer, self.nl_atom_numbers, self.nl_atom_serial, + self.uint_dr_to_dr_cof, self.excluded_list_start, self.excluded_list, + self.excluded_numbers, self.need_refresh_flag, self.refresh_count) # Done + self.LJ_Update_Volume() + self.PME_Update_Volume(factor) + self.Simple_Constrain_Update_Volume() + # self.mol_map.Update_Volume(self.md_info.sys.box_length) + + def Main_Volume_Change_Largely(self): + """Main_Volume_Change_Largely""" + # re-initialize neighbor_list and pme + _ = self.neighbor_list_update_forced_update(self.atom_numbers_in_grid_bucket, self.bucket, + self.crd, self.box_length, self.grid_N, + self.grid_length_inverse, self.atom_in_grid_serial, + self.old_crd, self.crd_to_uint_crd_cof, self.uint_crd, + self.pointer, self.nl_atom_numbers, self.nl_atom_serial, + self.uint_dr_to_dr_cof, self.excluded_list_start, + self.excluded_list, + self.excluded_numbers, self.need_refresh_flag, + self.refresh_count) + + def Check_MC_Barostat_Accept(self): + """Check_MC_Barostat_Accept""" + self.total_count = self.total_count + 1 + rand_num = self.random((1, 1)) + if rand_num[0] < self.mc_baro_accept_possibility: + self.reject = 0 + self.accept_count += 1 + else: + self.reject = 1 + return self.reject + + def Delta_V_Max_Update(self): + """Delta_V_Max_Update""" + if self.total_count % self.check_interval == 0: + self.accept_rate = 100.0 * self.accept_count / self.total_count + if self.accept_rate < self.accept_rate_low: + self.total_count = 0 + self.accept_count = 0 + self.DeltaV_max = self.DeltaV_max * 0.9 + if self.accept_rate > self.accept_rate_high: + self.total_count = 0 + self.accept_count = 0 + self.DeltaV_max = self.DeltaV_max * 1.1 + + def Main_iteration_presssure(self, steps, force): + """Main_iteration_presssure""" + if self.mc_baro_is_initialized and steps % self.mc_baro.update_interval == 0: + # old energy + self.mc_baro_energy_old = self.potential + self.frc_backup = self.frc + self.crd_backup = self.crd + self.Volume_Change_Attempt(self.box_length, 200) + + # change coordinates + if self.is_molecule_map_output: + nowrap_crd = self.Calculate_No_Wrap_Crd() + self.crd, _ = self.Residue_Crd_Map(nowrap_crd) + _ = self.refresh_boxmaptimes(self.crd, self.old_crd, 1.0 / self.box_length, self.box_map_times) + else: + self.crd = self.crd * self.crd_scale_factor # scale list + + # change volume + self.Main_Volume_Change(self.crd_scale_factor) + self.system_reinitializing_count += 1 + + # new energy + _ = self.Simulation_Caculate_Force(self.uint_crd, self.uint_dr_to_dr_cof, self.nl_atom_numbers, + self.nl_atom_serial) + + self.energy_new = self.potential + + # calculate accepted rate + if self.scale_coordinate_by_residue: + self.extra_term = self.target_pressure * self.DeltaV - \ + self.residue_numbers * self.CONSTANT_kB * \ + self.target_temperature * self.log(self.VDevided) + else: + self.extra_term = self.target_pressure * self.DeltaV - \ + self.atom_numbers * self.CONSTANT_kB * \ + self.target_temperature * self.log(self.VDevided) + + self.mc_baro_accept_possibility = self.mc_baro_energy_new - self.mc_baro_energy_old + self.extra_term + self.mc_baro.mc_baro_accept_possibility = self.exp( + -self.mc_baro_accept_possibility / (self.CONSTANT_kB * self.target_temperature)) + + # check if accepted + if self.Check_MC_Barostat_Accept(): + # if accept, refresh + self.crd_scale_factor = 1.0 / self.crd_scale_factor + self.crd = self.crd_backup + self.Main_Volume_Change(self.crd_scale_factor) + self.system_reinitializing_count += 1 + _ = self.neighbor_list_update_mc(self.atom_numbers_in_grid_bucket, self.bucket, + self.crd, self.box_length, self.grid_N, + self.grid_length_inverse, self.atom_in_grid_serial, + self.old_crd, self.crd_to_uint_crd_cof, self.uint_crd, + self.pointer, self.nl_atom_numbers, self.nl_atom_serial, + self.uint_dr_to_dr_cof, self.excluded_list_start, self.excluded_list, + self.excluded_numbers, self.need_refresh_flag, + self.refresh_count) + self.frc = force + self.frc = self.frc_backup + + # reinitialized + if self.system_reinitializing_count >= 20000 or (not self.reject and ( + self.mc_baro_newV > 1.331 * self.mc_baro_V0 or self.mc_baro_newV < 0.729 * self.mc_baro.V0)): + self.Main_Volume_Change_Largely() + self.mc_baro_V0 = self.mc_baro_newV + self.system_reinitializing_count = self.zero_fp_tensor + self.Delta_V_Max_Update() + + def Constrain(self): + """Constrain""" + constrain_frc = self.zero_frc + for _ in range(self.iteration_numbers): + test_uint_crd = self.refresh_uint_crd(self.crd, self.quarter_crd_to_uint_crd_cof, constrain_frc, + self.mass_inverse) + if self.need_pressure: + force, _ = self.constrain_force_cycle_with_virial(test_uint_crd, self.uint_dr_to_dr_cof, + self.last_pair_dr, self.atom_i_serials, + self.atom_j_serials, self.constant_rs, + self.constrain_ks) + else: + force = self.constrain_force_cycle(test_uint_crd, self.uint_dr_to_dr_cof, self.last_pair_dr, + self.atom_i_serials, + self.atom_j_serials, self.constant_rs, self.constrain_ks) + constrain_frc = constrain_frc + force + + res = self.refresh_crd_vel(self.crd, self.velocity, constrain_frc, self.mass_inverse) + crd = self.depend(self.crd, res) + vel = self.depend(self.velocity, res) + + return crd, vel, res + + def Main_Iteration(self, steps, force): + '''Main_Iteration''' + # self.Main_iteration_presssure(steps, force) + # Remember_Last_Coordinates + # pressure control 1 + if self.simple_constrain_is_initialized: + self.last_pair_dr = self.last_crd_to_dr(self.crd, self.quarter_crd_to_uint_crd_cof, self.uint_dr_to_dr_cof, + self.atom_i_serials, + self.atom_j_serials, self.constant_rs, self.constrain_ks) + + if self.mode == 0: # NVE + self.velocity, self.crd, _ = self.Simulation_MDIterationLeapFrog(force) + elif self.mode == -1: # Minimization + _ = self.Simulation_MDInformationGradientDescent(force) + else: + if self.liujian_info_is_initialized: + self.velocity, self.crd, _ = self.Simulation_MDIterationLeapFrog_Liujian(self.mass_inverse, + self.sqrt_mass, self.crd, + force, + self.rand_state, + self.random_force) + + if self.simple_constrain_is_initialized: + self.crd, self.velocity, res1 = self.Constrain() + else: + res1 = self.zero_fp_tensor + + # MD_Information_Crd_To_Uint_Crd + self.uint_crd = self.crd_to_uint_crd_quarter(self.quarter_crd_to_uint_crd_cof, self.crd) + res2 = self.neighbor_list_update(self.atom_numbers_in_grid_bucket, + self.bucket, + self.crd, + self.box_length, + self.grid_N, + self.grid_length_inverse, + self.atom_in_grid_serial, + self.old_crd, + self.crd_to_uint_crd_cof, + self.uint_crd, + self.pointer, + self.nl_atom_numbers, + self.nl_atom_serial, + self.uint_dr_to_dr_cof, + self.excluded_list_start, + self.excluded_list, + self.excluded_numbers, + self.need_refresh_flag, + self.refresh_count) + + res3 = self.refresh_boxmaptimes(self.crd, self.old_crd, 1.0 / self.box_length, self.box_map_times) + + return self.velocity, self.crd, res1, res2, res3 + + def Calculate_No_Wrap_Crd(self): + """Calculate_No_Wrap_Crd""" + nowrap_crd = self.box_map_times * self.box_length + self.crd + return nowrap_crd + + def Residue_Crd_Map(self, nowrap_crd): + """Residue_Crd_Map""" + center_of_mass = self.getcenterofmass(self.res_start, self.res_end, nowrap_crd, self.mass, + self.res_mass_inverse) + + res = self.mapcenterofmass(self.res_start, self.res_end, center_of_mass, self.box_length, nowrap_crd, self.crd) + + return self.crd, res + def construct(self, step, print_step): '''construct''' - self.last_crd = self.crd - res = self.neighbor_list_update(self.atom_numbers_in_grid_bucket, - self.bucket, - self.crd, - self.box_length, - self.grid_N, - self.grid_length_inverse, - self.atom_in_grid_serial, - self.old_crd, - self.crd_to_uint_crd_cof, - self.uint_crd, - self.pointer, - self.nl_atom_numbers, - self.nl_atom_serial, - self.uint_dr_to_dr_cof, - self.excluded_list_start, - self.excluded_list, - self.excluded_numbers, - self.need_refresh_flag, - self.refresh_count) - uint_crd = self.Simulation_Beforce_Caculate_Force() - force = self.Simulation_Caculate_Force(uint_crd, self.uint_dr_to_dr_cof, self.nl_atom_numbers, + # self.last_crd = self.crd + if step == 0: + res = self.neighbor_list_update_forced_update(self.atom_numbers_in_grid_bucket, + self.bucket, + self.crd, + self.box_length, + self.grid_N, + self.grid_length_inverse, + self.atom_in_grid_serial, + self.old_crd, + self.crd_to_uint_crd_cof, + self.uint_crd, + self.pointer, + self.nl_atom_numbers, + self.nl_atom_serial, + self.uint_dr_to_dr_cof, + self.excluded_list_start, + self.excluded_list, + self.excluded_numbers, + self.need_refresh_flag, + self.refresh_count) + else: + res = self.zero_fp_tensor + force = self.Simulation_Caculate_Force(self.uint_crd, self.uint_dr_to_dr_cof, self.nl_atom_numbers, self.nl_atom_serial) + if step == 0: + self.rand_state = self.setup_random_state() + + self.velocity, self.crd, res1, res2, res3 = self.Main_Iteration(step + 1, force) + temperature = self.Simulation_Temperature() if print_step == 0: bond_energy_sum, angle_energy_sum, dihedral_energy_sum, nb14_lj_energy_sum, nb14_cf_energy_sum, \ - lj_energy_sum, ee_ene, total_energy = self.Simulation_Caculate_Energy(uint_crd, self.uint_dr_to_dr_cof) + lj_energy_sum, ee_ene, total_energy = self.Simulation_Caculate_Energy(self.uint_crd, self.uint_dr_to_dr_cof) else: bond_energy_sum = self.zero_fp_tensor angle_energy_sum = self.zero_fp_tensor @@ -442,12 +960,5 @@ class Simulation(nn.Cell): lj_energy_sum = self.zero_fp_tensor ee_ene = self.zero_fp_tensor total_energy = self.zero_fp_tensor - temperature = self.Simulation_Temperature() - if step == 0: - self.rand_state = self.setup_random_state() - self.velocity, self.crd, _ = self.Simulation_MDIterationLeapFrog_Liujian(self.mass_inverse, - self.sqrt_mass, self.crd, force, - self.rand_state, - self.random_force) return temperature, total_energy, bond_energy_sum, angle_energy_sum, dihedral_energy_sum, nb14_lj_energy_sum, \ - nb14_cf_energy_sum, lj_energy_sum, ee_ene, res + nb14_cf_energy_sum, lj_energy_sum, ee_ene, res, res1, res2, res3 diff --git a/model_zoo/research/hpc/sponge/src/system_information.py b/model_zoo/research/hpc/sponge/src/system_information.py new file mode 100644 index 00000000000..ad99526c824 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/system_information.py @@ -0,0 +1,290 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''System information''' +import numpy as np + + +class periodic_box_condition_information: + """periodic_box_condition_information""" + + def __init__(self, controller, box_length): + CONSTANT_UINT_MAX_FLOAT = 4294967296.0 + self.crd_to_uint_crd_cof = np.array([CONSTANT_UINT_MAX_FLOAT / box_length[0], + CONSTANT_UINT_MAX_FLOAT / box_length[1], + CONSTANT_UINT_MAX_FLOAT / box_length[2]]) + self.quarter_crd_to_uint_crd_cof = 0.25 * self.crd_to_uint_crd_cof + self.uint_dr_to_dr_cof = 1.0 / self.crd_to_uint_crd_cof + + +class system_information: + """system_information""" + + def __init__(self, controller, md_info): + CONSTANT_PRES_CONVERTION_INVERSE = 0.00001439506089041446 + self.md_info = md_info + self.box_length = self.md_info.box_length + self.steps = 0 + self.step_limit = 1000 if "step_limit" not in controller.Command_Set else int( + controller.Command_Set["step_limit"]) + self.target_temperature = 300.0 if "target_temperature" not in controller.Command_Set else float( + controller.Command_Set["target_temperature"]) + if md_info.mode == 2 and "target_pressure" in controller.Command_Set: + self.target_pressure = float(controller.Command_Set["target_pressure"]) + else: + self.target_pressure = 1 + self.target_pressure *= CONSTANT_PRES_CONVERTION_INVERSE + self.d_virial = 0 + self.d_pressure = 0 + self.d_temperature = 0 + self.d_potential = 0 + self.d_sum_of_atom_ek = 0 + self.freedom = 3 * md_info.atom_numbers + + +class non_bond_information: + """system_information""" + + def __init__(self, controller, md_info): + self.md_info = md_info + self.skin = 2.0 if "skin" not in controller.Command_Set else float(controller.Command_Set["skin"]) + print(" skin set to %.2f Angstram" % (self.skin)) + self.cutoff = 10.0 if "cutoff" not in controller.Command_Set else float(controller.Command_Set["cutoff"]) + self.atom_numbers = self.md_info.atom_numbers + self.excluded_atom_numbers = 0 + self.h_excluded_list_start = [] + self.h_excluded_numbers = [] + self.h_excluded_list = [] + if controller.amber_parm is not None: + file_path = controller.amber_parm + self.read_information_from_amberfile(file_path) + else: + self.read_exclude_file(controller) + + def read_exclude_file(self, controller): + """read_exclude_file""" + if "exclude_in_file" in controller.Command_Set: + print(" Start reading excluded list:") + path = controller.Command_Set["exclude_in_file"] + file = open(path, 'r') + context = file.readlines() + atom_numbers, self.excluded_atom_numbers = list(map(int, context[0].strip().split())) + if self.md_info.atom_numbers > 0 and (atom_numbers != self.md_info.atom_numbers): + print(" Error: atom_numbers is not equal: ", atom_numbers, self.md_info.atom_numbers) + exit(1) + else: + self.md_info.atom_numbers = atom_numbers + count = 0 + for idx, val in enumerate(context): + if idx > 0: + el = list(map(int, val.strip().split())) + if el[0] == 1 and -1 in el: + self.h_excluded_numbers.append(0) + else: + self.h_excluded_numbers.append(el[0]) + self.h_excluded_list_start.append(count) + if el: + self.h_excluded_list.extend(el[1:]) + count += el[0] + print(" End reading excluded list") + file.close() + else: + print(" Set all atom exclude no atoms as default") + count = 0 + for i in range(self.md_info.atom_numbers): + self.h_excluded_numbers[i] = 0 + self.h_excluded_list_start[i] = count + for _ in range(self.h_excluded_numbers[i]): + self.h_excluded_list[count] = 0 + count += 1 + print(" End reading charge") + + def read_information_from_amberfile(self, file_path): + '''read amber file''' + file = open(file_path, 'r') + context = file.readlines() + file.close() + self.h_excluded_list_start = [0] * self.atom_numbers + self.h_excluded_numbers = [0] * self.atom_numbers + + for idx, val in enumerate(context): + if idx < len(context) - 1: + if "%FLAG POINTERS" in val + context[idx + 1] and "%FORMAT(10I8)" in val + context[idx + 1]: + start_idx = idx + 2 + count = 0 + value = list(map(int, context[start_idx].strip().split())) + information = [] + information.extend(value) + while count < 11: + start_idx += 1 + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + self.excluded_atom_numbers = information[10] + print("excluded atom numbers ", self.excluded_atom_numbers) + break + for idx, val in enumerate(context): + if "%FLAG NUMBER_EXCLUDED_ATOMS" in val: + count = 0 + start_idx = idx + information = [] + while count < self.atom_numbers: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + count = 0 + for i in range(self.atom_numbers): + self.h_excluded_numbers[i] = information[i] + self.h_excluded_list_start[i] = count + count += information[i] + break + + total_count = sum(self.h_excluded_numbers) + self.h_excluded_list = [] + for idx, val in enumerate(context): + if "%FLAG EXCLUDED_ATOMS_LIST" in val: + count = 0 + start_idx = idx + information = [] + while count < total_count: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + + count = 0 + for i in range(self.atom_numbers): + tmp_list = [] + if self.h_excluded_numbers[i] == 1: + tmp_list.append(information[count] - 1) + if information[count] == 0: + self.h_excluded_numbers[i] = 0 + count += 1 + else: + for _ in range(self.h_excluded_numbers[i]): + tmp_list.append(information[count] - 1) + + count += 1 + tmp_list = sorted(tmp_list) + self.h_excluded_list.extend(tmp_list) + break + + +class NVE_iteration: + """NVE_iteration""" + + def __init__(self, controller, md_info): + self.max_velocity = -1 if "nve_velocity_max" not in controller.Command_Set else float( + controller.Command_Set["nve_velocity_max"]) + + +class residue_information: + """residue_information""" + + def __init__(self, controller, md_info): + self.md_info = md_info + self.residue_numbers = 0 + self.h_mass = [] + self.h_mass_inverse = [] + self.h_res_start = [] + self.h_res_end = [] + self.momentum = [] + self.center_of_mass = [] + self.sigma_of_res_ek = 0 + self.res_ek_energy = 0 + self.sigma_of_res_ek = 0 + self.is_initialized = 0 + print(" Start reading residue list:") + if "residue_in_file" in controller.Command_Set: + self.read_residule_file(controller) + elif "amber_parm7" in controller.Command_Set: + self.residue_numbers = self.md_info.residue_numbers + self.h_res_start = md_info.h_res_start + self.h_res_end = md_info.h_res_end + self.is_initialized = 1 + self.read_res_mass() + else: + self.residue_numbers = md_info.atom_numbers + self.h_res_start = list(range(self.residue_numbers)) + self.h_res_end = list(range(1, self.residue_numbers + 1)) + self.is_initialized = 1 + print(" End reading residue list") + + def read_res_mass(self): + """ Read_AMBER_Parm7 """ + if self.md_info.h_mass: + for i in range(self.residue_numbers): + temp_mass = 0 + for j in range(self.h_res_start[i], self.h_res_end[i]): + temp_mass += self.md_info.h_mass[j] + self.h_mass.append(temp_mass) + if temp_mass == 0: + self.h_mass_inverse.append(0) + else: + self.h_mass_inverse.append(1.0 / temp_mass) + else: + print(" Error: atom mass should be initialized before residue mass") + exit(1) + + def read_residule_file(self, controller): + """read_residule_file""" + if "residue_in_file" in controller.Command_Set: + path = controller.Command_Set["residue_in_file"] + file = open(path, 'r') + context = file.readlines() + atom_numbers, self.residue_numbers = list(map(int, context[0].strip().split())) + print(" residue_numbers is ", self.residue_numbers) + # self.md_info.residue_numbers = self.residue_numbers + if self.md_info.atom_numbers > 0 and (atom_numbers != self.md_info.atom_numbers): + print(" Error: atom_numbers is not equal: ", atom_numbers, self.md_info.atom_numbers) + exit(1) + else: + self.md_info.atom_numbers = atom_numbers + print(" residue_numbers is ", self.residue_numbers) + + count = 0 + for idx, val in enumerate(context): + if idx > 0: + self.h_res_start.append(count) + temp = int(val.strip()) + count += temp + self.h_res_end.append(count) + print(" End reading excluded list") + file.close() + self.is_initialized = 1 + if self.is_initialized: + self.read_res_mass() + + +class trajectory_output: + """trajectory_output""" + + def __init__(self, controller, md_info): + self.current_crd_synchronized_step = 0 + self.is_molecule_map_output = 0 + if "molecule_map_output" in controller.Command_Set: + self.is_molecule_map_output = int(controller.Command_Set["molecule_map_output"]) + self.amber_irest = -1 + self.write_trajectory_interval = 1000 if "write_information_interval" not in controller.Command_Set else int( + controller.Command_Set["write_information_interval"]) + self.write_restart_file_interval = self.write_trajectory_interval if "write_restart_file_interval" not in \ + controller.Command_Set else \ + int(controller.Command_Set["write_restart_file_interval"]) diff --git a/model_zoo/research/hpc/sponge/src/vatom.py b/model_zoo/research/hpc/sponge/src/vatom.py new file mode 100644 index 00000000000..05b90a7e174 --- /dev/null +++ b/model_zoo/research/hpc/sponge/src/vatom.py @@ -0,0 +1,287 @@ +# Copyright 2021 Huawei Technologies Co., Ltd +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ============================================================================ +'''virtual_information''' + + +class VIRTUAL_TYPE_INFROMATION: + """VIRTUAL_LAYER_INFORMATION""" + + def __init__(self): + self.virtual_numbers = 0 + self.virtual_type = [] + + +class VIRTUAL_LAYER_INFORMATION: + '''VIRTUAL_LAYER_INFORMATION''' + + def __init__(self): + self.v0_info = VIRTUAL_TYPE_INFROMATION() + self.v1_info = VIRTUAL_TYPE_INFROMATION() + self.v2_info = VIRTUAL_TYPE_INFROMATION() + self.v3_info = VIRTUAL_TYPE_INFROMATION() + + +class Virtual_Information: + '''virtual_information''' + + def __init__(self, controller, md_info, system_freedom): + self.module_name = "virtual_atom" + self.atom_numbers = md_info.atom_numbers + self.virtual_level = [0] * self.atom_numbers + self.virtual_layer_info = [] + self.system_freedom = system_freedom + self.max_level = 1 + self.is_initialized = 0 + name = self.module_name + "_in_file" + if name in controller.Command_Set: + self.virtual_layer_info_2 = [] + path = controller.Command_Set[name] + print(" Start reading virtual levels\n") + self.read_in_file(path) + self.max_level, self.total_virtual_atoms = self.level_init() + self.system_freedom -= 3 * self.total_virtual_atoms + print(" Virtual Atoms Max Level is ", self.max_level) + print(" Virtual Atoms Number is ", self.total_virtual_atoms) + print(" End reading virtual levels") + self.read_in_file_second(path) + self.read_in_file_third(path) + self.is_initialized = 1 + + if controller.amber_parm is not None: + file_path = controller.amber_parm + self.read_information_from_amberfile(file_path) + + def read_in_file(self, path): + """read_in_file""" + file = open(path, 'r') + context = file.readlines() + line_numbers = 0 + for _, val in enumerate(context): + line_numbers += 1 + tl = list(map(float, val.strip().split())) + virtual_type, virtual_atom = int(tl[0]), int(tl[1]) + if virtual_type == 0: + temp, _ = [int(tl[2])], [tl[3]] + self.virtual_level[virtual_atom] = self.virtual_level[temp[0]] + 1 + if virtual_type == 1: + temp, _ = [int(tl[2]), int(tl[3])], [tl[4]] + self.virtual_level[virtual_atom] = max(self.virtual_level[temp[0]], self.virtual_level[temp[1]]) + 1 + if virtual_type == 2: + temp, _ = [int(tl[2]), int(tl[3]), int(tl[4])], [tl[5], tl[6]] + self.virtual_level[virtual_atom] = max(self.virtual_level[temp[0]], + self.virtual_level[temp[1]], + self.virtual_level[temp[2]]) + 1 + if virtual_type == 3: + temp, _ = [int(tl[2]), int(tl[3]), int(tl[4])], [tl[5], tl[6]] + self.virtual_level[virtual_atom] = max(self.virtual_level[temp[0]], + self.virtual_level[temp[1]], + self.virtual_level[temp[2]]) + 1 + if virtual_type > 3 or virtual_type < 0: + print(" Error: can not parse line #{} because {} is not a proper type for virtual atoms.".format( + line_numbers, virtual_type)) + exit(1) + file.close() + + def level_init(self): + """level_init""" + max_level = 0 + total_virtual_atoms = 0 + for i in range(self.atom_numbers): + vli = self.virtual_level[i] + if vli > 0: + total_virtual_atoms += 1 + if vli > max_level: + for _ in range(vli - max_level): + virtual_layer = VIRTUAL_LAYER_INFORMATION() + # v0_info.virtual_numbers, v1_info.virtual_numbers + # v2_info.virtual_numbers, v3_info.virtual_numbers + self.virtual_layer_info.append(virtual_layer) + max_level = vli + return max_level, total_virtual_atoms + + def read_in_file_second(self, path): + """read_in_file_second""" + print(" Start reading virtual type numbers in different levels") + file = open(path, 'r') + context = file.readlines() + line_numbers = 0 + for _, val in enumerate(context): + line_numbers += 1 + tl = list(map(float, val.strip().split())) + virtual_type, virtual_atom = int(tl[0]), int(tl[1]) + temp_vl = self.virtual_layer_info[self.virtual_level[virtual_atom] - 1] + if virtual_type == 0: + temp_vl.v0_info.virtual_numbers += 1 + if virtual_type == 1: + temp_vl.v1_info.virtual_numbers += 1 + if virtual_type == 2: + temp_vl.v2_info.virtual_numbers += 1 + if virtual_type == 3: + temp_vl.v3_info.virtual_numbers += 1 + + self.virtual_layer_info[self.virtual_level[virtual_atom] - 1] = temp_vl + print(" End reading virtual type numbers in different levels") + + def read_in_file_third(self, path): + """read_in_file_third""" + print(" Start reading information for every virtual atom") + file = open(path, 'r') + context = file.readlines() + line_numbers = 0 + count0, count1, count2, count3 = 0, 0, 0, 0 + temp_v = VIRTUAL_LAYER_INFORMATION() + self.virtual_layer_info_2 = [0] * len(self.virtual_layer_info) + for _, val in enumerate(context): + line_numbers += 1 + tl = list(map(float, val.strip().split())) + virtual_type, virtual_atom = int(tl[0]), int(tl[1]) + temp_vl = self.virtual_layer_info[self.virtual_level[virtual_atom] - 1] + if virtual_type == 0: + temp_vl.v0_info.virtual_type.append([int(tl[1]), int(tl[2]), 2 * tl[3]]) + # virtual_atom, from_1, from_2, from_3, h_double + count0 += 1 + if virtual_type == 1: + temp_v.v1_info.virtual_type.append([int(tl[1]), int(tl[2]), int(tl[3]), tl[4]]) + # virtual_atom, from_1, from_2, from_3, a + count1 += 1 + if virtual_type == 2: + temp_v.v2_info.virtual_type.append([int(tl[1]), int(tl[2]), int(tl[3]), int(tl[4]), tl[5], tl[6]]) + # virtual_atom, from_1, from_2, from_3, a, b + count2 += 1 + if virtual_type == 3: + temp_v.v3_info.virtual_type.append([int(tl[1]), int(tl[2]), int(tl[3]), int(tl[4]), tl[5], tl[6]]) + # virtual_atom, from_1, from_2, from_3, d, k + count3 += 1 + self.virtual_layer_info_2[self.virtual_level[virtual_atom] - 1] = temp_v + file.close() + print(" End reading information for every virtual atom") + + def read_information_from_amberfile(self, file_path): + '''read amber file''' + file = open(file_path, 'r') + context = file.readlines() + file.close() + for idx, val in enumerate(context): + if idx < len(context) - 1: + if "%FLAG POINTERS" in val + context[idx + 1] and "%FORMAT(10I8)" in val + context[idx + 1]: + start_idx = idx + 2 + count = 0 + value = list(map(int, context[start_idx].strip().split())) + self.angle_with_H_numbers = value[4] + self.angle_without_H_numbers = value[5] + self.angle_numbers = self.angle_with_H_numbers + self.angle_without_H_numbers + information = [] + information.extend(value) + while count < 15: + start_idx += 1 + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + self.angle_type_numbers = information[16] + print("angle type numbers ", self.angle_type_numbers) + break + + self.h_atom_a = [0] * self.angle_numbers + self.h_atom_b = [0] * self.angle_numbers + self.h_atom_c = [0] * self.angle_numbers + self.h_type = [0] * self.angle_numbers + angle_count = 0 + for idx, val in enumerate(context): + if "%FLAG ANGLES_INC_HYDROGEN" in val: + count = 0 + start_idx = idx + information = [] + while count < 4 * self.angle_with_H_numbers: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + for _ in range(self.angle_with_H_numbers): + self.h_atom_a[angle_count] = information[angle_count * 4 + 0] / 3 + self.h_atom_b[angle_count] = information[angle_count * 4 + 1] / 3 + self.h_atom_c[angle_count] = information[angle_count * 4 + 2] / 3 + self.h_type[angle_count] = information[angle_count * 4 + 3] - 1 + angle_count += 1 + + break + + for idx, val in enumerate(context): + if "%FLAG ANGLES_WITHOUT_HYDROGEN" in val: + count = 0 + start_idx = idx + information = [] + while count < 4 * self.angle_without_H_numbers: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(int, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + for _ in range(self.angle_without_H_numbers): + self.h_atom_a[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 0] / 3 + self.h_atom_b[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 1] / 3 + self.h_atom_c[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 2] / 3 + self.h_type[angle_count] = information[(angle_count - self.angle_with_H_numbers) * 4 + 3] - 1 + angle_count += 1 + break + self.processor(context, angle_count) + + def processor(self, context, angle_count): + ''' processor ''' + self.type_k = [0] * self.angle_type_numbers + for idx, val in enumerate(context): + if "%FLAG ANGLE_FORCE_CONSTANT" in val: + count = 0 + start_idx = idx + information = [] + while count < self.angle_type_numbers: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(float, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + self.type_k = information[:self.angle_type_numbers] + break + + self.type_theta0 = [0] * self.angle_type_numbers + for idx, val in enumerate(context): + if "%FLAG ANGLE_EQUIL_VALUE" in val: + count = 0 + start_idx = idx + information = [] + while count < self.angle_type_numbers: + start_idx += 1 + if "%FORMAT" in context[start_idx]: + continue + else: + value = list(map(float, context[start_idx].strip().split())) + information.extend(value) + count += len(value) + self.type_theta0 = information[:self.angle_type_numbers] + break + if self.angle_numbers != angle_count: + print("angle count %d != angle_number %d ", angle_count, self.angle_numbers) + + self.h_angle_k = [] + self.h_angle_theta0 = [] + for i in range(self.angle_numbers): + self.h_angle_k.append(self.type_k[self.h_type[i]]) + self.h_angle_theta0.append(self.type_theta0[self.h_type[i]])