!21243 sponge case update

Merge pull request !21243 from jiahongQian/master_case

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()

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_a[angle_count] = int(information[angle_count * 4 + 0] / 3)
-                    self.h_atom_b[angle_count] = information[angle_count * 4 + 1] / 3
+                    self.h_atom_b[angle_count] = int(information[angle_count * 4 + 1] / 3)
-                    self.h_atom_c[angle_count] = information[angle_count * 4 + 2] / 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_a[angle_count] = int(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_b[angle_count] = int(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_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

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")

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_a[i] = int(information[3 * i + 0] / 3)
-                    self.h_atom_b[i] = information[3 * i + 1] / 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_a[i] = int(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_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]

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]

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.h_pk = []
-        self.gamc = []
+        self.h_gamc = []
-        self.gams = []
+        self.h_gams = []
-        self.pn = []
+        self.h_pn = []
-        self.ipn = []
+        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.h_gams.append(tmpf2 * self.h_pk[i])
-                    self.pn.append(abs(self.pn_type[tmpi]))
+                    self.h_pn.append(abs(self.pn_type[tmpi]))
-                    self.ipn.append(int(self.pn[i] + 0.001))
+                    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.h_gams.append(tmpf2 * self.h_pk[i])
-                    self.pn.append(abs(self.pn_type[tmpi]))
+                    self.h_pn.append(abs(self.pn_type[tmpi]))
-                    self.ipn.append(int(self.pn[i] + 0.001))
+                    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:

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(
+        self.gamma_ln = 1.0
-            controller.Command_Set["langevin_gamma"])
+        if "gamma" in controller.Command_Set:
-        self.rand_seed = 1 if "langevin_seed" not in controller.Command_Set else float(
+            self.gamma_ln = float(controller.Command_Set["gamma"])
-            controller.Command_Set["langevin_seed"])
+        if "langevin_liu_gamma" in controller.Command_Set:
-        self.max_velocity = 10000.0 if "velocity_max" not in controller.Command_Set else float(
+            self.gamma_ln = float(controller.Command_Set["langevin_liu_gamma"])
-            controller.Command_Set["velocity_max"])
+        print("    langevin_liu_gamma is ", self.gamma_ln)
-        assert self.max_velocity > 0
+
-        print("target temperature is ", self.target_temperature)
+        self.random_seed = 1 if "seed" not in controller.Command_Set else int(
-        print("friction coefficient is ", self.gamma_ln, "ps^-1")
+            controller.Command_Set["seed"])
-        print("random seed is ", self.rand_seed)
+
-        self.dt = float(controller.Command_Set["dt"])
+        print("    target temperature is {} K".format(self.target_temperature))
-        self.dt *= self.CONSTANT_TIME_CONVERTION
+        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'''

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

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")

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],
+        self.sys = system_information(controller, self)
-                                    CONSTANT_UINT_MAX_FLOAT / self.box_length[1],
+        self.nb = non_bond_information(controller, self)
-                                    CONSTANT_UINT_MAX_FLOAT / self.box_length[2]]
+        self.output = trajectory_output(controller, self)
-        self.uint_dr_to_dr_cof = [1.0 / self.crd_to_uint_crd_cof[0], 1.0 / self.crd_to_uint_crd_cof[1],
+        self.nve = NVE_iteration(controller, self)
-                                  1.0 / self.crd_to_uint_crd_cof[2]]
+        self.res = residue_information(controller, self)
-        self.density *= 1e24 / 6.023e23 / (self.box_length[0] * self.box_length[1] * self.box_length[2])
+        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])

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_a = self.h_atom_a[:self.nb14_numbers]
-        self.h_atom_b = self.h_atom_b[: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_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_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'''

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(
+        self.CONSTANT_UINT_MAX_FLOAT = 4294967296.0
-            controller.Command_Set["neighbor_list_refresh_interval"])
+        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]]

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.module_name = "PME"
-        self.tolerance = 0.00001 if "PME_Direct_Tolerance" not in controller.Command_Set else float(
+        self.CONSTANT_Pi = 3.1415926535897932
-            controller.Command_Set["PME_Direct_Tolerance"])
+        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'''

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())))

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()

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
+import numpy as np
-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
 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.crd_to_uint_crd_cof = Tensor(np.asarray(self.md_info.pbc.crd_to_uint_crd_cof, np.float32), mstype.float32)
-        self.uint_dr_to_dr_cof = Parameter(
+        self.quarter_crd_to_uint_crd_cof = Tensor(np.asarray(self.md_info.pbc.quarter_crd_to_uint_crd_cof, np.float32),
-            Tensor(np.asarray(self.md_info.uint_dr_to_dr_cof, np.float32), mstype.float32), requires_grad=False)
+                                                  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),
+        self.mass = Tensor(self.md_info.h_mass, mstype.float32)
-                                  requires_grad=False)
         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.pk = Tensor(np.asarray(self.dihedral.h_pk, np.float32), mstype.float32)
-        self.gamc = Tensor(np.asarray(self.dihedral.gamc, np.float32), mstype.float32)
+        self.gamc = Tensor(np.asarray(self.dihedral.h_gamc, np.float32), mstype.float32)
-        self.gams = Tensor(np.asarray(self.dihedral.gams, np.float32), mstype.float32)
+        self.gams = Tensor(np.asarray(self.dihedral.h_gams, np.float32), mstype.float32)
-        self.pn = Tensor(np.asarray(self.dihedral.pn, np.float32), mstype.float32)
+        self.pn = Tensor(np.asarray(self.dihedral.h_pn, np.float32), mstype.float32)
-        self.ipn = Tensor(np.asarray(self.dihedral.ipn, np.int32), mstype.int32)
+        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_start = Tensor(np.asarray(self.md_info.nb.h_excluded_list_start, np.int32), mstype.int32)
-        self.excluded_list = Tensor(np.asarray(self.nb_info.excluded_list, 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.nb_info.excluded_numbers, 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.atom_LJ_type = Tensor(self.LJ_info.atom_LJ_type, mstype.int32)
-        self.LJ_A = Tensor(np.asarray(self.LJ_info.LJ_A, dtype=np.float32), mstype.float32)
+        self.LJ_A = Tensor(self.LJ_info.h_LJ_A, mstype.float32)
-        self.LJ_B = Tensor(np.asarray(self.LJ_info.LJ_B, dtype=np.float32), 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,
+        self.md_iteration_leap_frog_liujian_with_max_vel = P.MDIterationLeapFrogLiujianWithMaxVel(self.atom_numbers,
-                                                              atom_numbers=self.atom_numbers, not_first_time=0,
+                                                                                                  self.half_dt, self.dt,
-                                                              nxy=self.nxy,
+                                                                                                  self.exp_gamma,
-                                                              excluded_atom_numbers=self.excluded_atom_numbers,
+                                                                                                  self.max_velocity)
-                                                              cutoff_square=self.cutoff_square,
+        self.neighbor_list_update = \
-                                                              half_skin_square=self.half_skin_square,
+            P.NeighborListUpdate(grid_numbers=self.grid_numbers,
-                                                              cutoff_with_skin=self.cutoff_with_skin,
+                                 atom_numbers=self.atom_numbers,
-                                                              half_cutoff_with_skin=self.half_cutoff_with_skin,
+                                 not_first_time=1, nxy=self.nxy,
-                                                              cutoff_with_skin_square=self.cutoff_with_skin_square,
+                                 excluded_atom_numbers=self.excluded_atom_numbers,
-                                                              refresh_interval=self.refresh_interval,
+                                 cutoff_square=self.cutoff_square,
-                                                              cutoff=self.cutoff, skin=self.skin,
+                                 half_skin_square=self.half_skin_square,
-                                                              max_atom_in_grid_numbers=self.max_atom_in_grid_numbers,
+                                 cutoff_with_skin=self.cutoff_with_skin,
-                                                              max_neighbor_numbers=self.max_neighbor_numbers)
+                                 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
+        self.uint_crd = self.crd_to_uint_crd_quarter(self.quarter_crd_to_uint_crd_cof, self.crd)
-        uint_crd = self.crd_to_uint_crd(crd_to_uint_crd_cof, self.crd)
+        return self.uint_crd
-        return 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,
+        uint_crd = self.Simulation_Beforce_Caculate_Force()
-                                                         self.bond_atom_b, self.bond_k, self.bond_r0)
+        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,
+        if self.pme_is_initialized:
-                                                           self.angle_atom_b, self.angle_atom_c,
+            pme_excluded_force = self.pme_excluded_force(uint_crd, scaler, self.charge, self.excluded_list_start,
-                                                           self.angle_k, self.angle_theta0)
+                                                         self.excluded_list, self.excluded_numbers)
 
-        dihedral_force, _ = self.dihedral_force_with_atom_energy(uint_crd, scaler,
+            pme_reciprocal_force = self.pme_reciprocal_force(uint_crd, self.charge)
-                                                                 self.dihedral_atom_a,
+            force = force + pme_excluded_force + pme_reciprocal_force
-                                                                 self.dihedral_atom_b,
+        if self.nb14_is_initialized:
-                                                                 self.dihedral_atom_c,
+            nb14_force, _ = self.nb14_force_with_atom_energy(uint_crd, self.atom_LJ_type, self.charge,
-                                                                 self.dihedral_atom_d, self.ipn,
+                                                             scaler, self.nb14_atom_a, self.nb14_atom_b,
-                                                                 self.pk, self.gamc, self.gams,
+                                                             self.lj_scale_factor, self.cf_scale_factor,
-                                                                 self.pn)
+                                                             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,
+        if self.bond_is_initialized:
-                                                         scaler, self.nb14_atom_a, self.nb14_atom_b,
+            bond_force, _, _ = self.bond_force_with_atom_energy_virial(uint_crd, scaler, self.bond_atom_a,
-                                                         self.lj_scale_factor, self.cf_scale_factor,
+                                                                       self.bond_atom_b, self.bond_k, self.bond_r0)
-                                                         self.LJ_A, self.LJ_B)
+            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,
+        if self.max_velocity <= 0:
-                                                  rand_state, random_frc)
+            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
+        # self.last_crd = self.crd
-        res = self.neighbor_list_update(self.atom_numbers_in_grid_bucket,
+        if step == 0:
-                                        self.bucket,
+            res = self.neighbor_list_update_forced_update(self.atom_numbers_in_grid_bucket,
-                                        self.crd,
+                                                          self.bucket,
-                                        self.box_length,
+                                                          self.crd,
-                                        self.grid_N,
+                                                          self.box_length,
-                                        self.grid_length_inverse,
+                                                          self.grid_N,
-                                        self.atom_in_grid_serial,
+                                                          self.grid_length_inverse,
-                                        self.old_crd,
+                                                          self.atom_in_grid_serial,
-                                        self.crd_to_uint_crd_cof,
+                                                          self.old_crd,
-                                        self.uint_crd,
+                                                          self.crd_to_uint_crd_cof,
-                                        self.pointer,
+                                                          self.uint_crd,
-                                        self.nl_atom_numbers,
+                                                          self.pointer,
-                                        self.nl_atom_serial,
+                                                          self.nl_atom_numbers,
-                                        self.uint_dr_to_dr_cof,
+                                                          self.nl_atom_serial,
-                                        self.excluded_list_start,
+                                                          self.uint_dr_to_dr_cof,
-                                        self.excluded_list,
+                                                          self.excluded_list_start,
-                                        self.excluded_numbers,
+                                                          self.excluded_list,
-                                        self.need_refresh_flag,
+                                                          self.excluded_numbers,
-                                        self.refresh_count)
+                                                          self.need_refresh_flag,
-        uint_crd = self.Simulation_Beforce_Caculate_Force()
+                                                          self.refresh_count)
-        force = self.Simulation_Caculate_Force(uint_crd, self.uint_dr_to_dr_cof, self.nl_atom_numbers,
+        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

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"])

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]])