1 0
import logging
2 0
import os
3

4 0
import numpy as np
5 0
from parmed.structure import Structure as ParmedStructureClass
6

7 0
from paprika.dummy import extract_dummy_atoms
8 0
from paprika.restraints.utils import get_bias_potential_type, parse_window
9 0
from paprika.utils import get_key, return_parmed_structure
10

11 0
logger = logging.getLogger(__name__)
12

13 0
_PI_ = np.pi
14

15

16 0
class Plumed:
17
    """
18
    This class converts restraints generated in pAPRika DAT_restraints into Plumed restraints.
19

20
    # TODO: possibly change this module to use the python wrapper of Plumed.
21

22
    Example:
23
    -------
24
        >>> plumed = Plumed()
25
        >>> plumed.file_name = 'plumed.dat'
26
        >>> plumed.path = './windows'
27
        >>> plumed.window_list = window_list
28
        >>> plumed.restraint_list = restraint_list
29
        >>> plumed.dump_to_file()
30

31
    plumed.dat:
32
        UNITS LENGTH=A ENERGY=kcal/mol TIME=ns
33
        # Collective variables
34
        c1 : DISTANCE ATOMS=175,150 NOPBC
35
        c2 : ANGLE    ATOMS=176,175,150 NOPBC
36
        c3 : ANGLE    ATOMS=175,150,165 NOPBC
37
        # Bias potential
38
        RESTRAINT   ARG=c7  AT=  6.000 KAPPA=  10.00
39
        RESTRAINT   ARG=c8  AT=  3.142 KAPPA= 200.00
40
        RESTRAINT   ARG=c9  AT=  3.142 KAPPA= 200.00
41
    """
42

43 0
    @property
44
    def file_name(self):
45
        """
46
        str: The Plumed file name where the restraints will be written to.
47
        """
48 0
        return self._file_name
49

50 0
    @file_name.setter
51 0
    def file_name(self, value: str):
52 0
        self._file_name = value
53

54 0
    @property
55
    def window_list(self):
56
        """
57
        list: The list of APR windows
58
        """
59 0
        return self._window_list
60

61 0
    @window_list.setter
62 0
    def window_list(self, value: list):
63 0
        self._window_list = value
64

65 0
    @property
66
    def restraint_list(self):
67
        """
68
        list: The list of restraints to be converted to Plumed.
69
        """
70 0
        return self._restraint_list
71

72 0
    @restraint_list.setter
73 0
    def restraint_list(self, value: list):
74 0
        self._restraint_list = value
75

76 0
    @property
77
    def path(self):
78
        """
79
        os.PathLike: The parent directory that contains the simulation windows.
80
        """
81 0
        return self._path
82

83 0
    @path.setter
84 0
    def path(self, value: str):
85 0
        self._path = value
86

87 0
    @property
88
    def uses_legacy_k(self):
89
        """
90
        bool: Option to specify whether the force constant parsed into DAT_restraint()
91
        is Amber-style or Gromacs/NAMD-style. Amber-style force constants have their value
92
        multiplied by a factor of 1/2 whereas Gromacs/NAMD-style does not. Plumed follows
93
        the Gromacs/NAMD-style for the force constant and the equations below demonstrates
94
        this point.
95

96
            * Amber:  U = K x (x - x0)²    * Plumed: U = 1/2 x k x (x - x0)²
97
            --> K(Amber) = 1/2 k(Plumed)
98

99
        i.e. "uses_legacy_k" is set to True (default) the force constants will be multiplied by 2.
100
        """
101 0
        return self._uses_legacy_k
102

103 0
    @uses_legacy_k.setter
104 0
    def uses_legacy_k(self, value: bool):
105 0
        self._uses_legacy_k = value
106

107 0
    @property
108
    def units(self):
109
        """
110
        dict: Dictionary of units for Plumed as strings. The dictionary requires the key values
111
        of 'energy', 'length, 'time'.
112
        """
113 0
        return self._units
114

115 0
    @units.setter
116 0
    def units(self, value: dict):
117 0
        self._units = value
118

119 0
    def __init__(self):
120 0
        self._file_name = "plumed.dat"
121 0
        self._restraint_list = None
122 0
        self._window_list = None
123 0
        self._path = "./"
124 0
        self._uses_legacy_k = True
125 0
        self.k_factor = 1.0
126 0
        self._units = None
127 0
        self.header_line = None
128 0
        self.group_index = None
129 0
        self.group_atoms = None
130

131 0
    def _initialize(self):
132
        # Set factor for spring constant
133 0
        if self.uses_legacy_k:
134 0
            self.k_factor = 2.0
135

136
        # Check units
137
        # NOTE: We have to resort to strings until (if) we migrate all quantities to Pint/SimTK units
138 0
        if self.units is None:
139 0
            self.units = {"energy": "kcal/mol", "length": "A", "time": "ns"}
140 0
        _check_plumed_units(self.units)
141

142
        # header line
143 0
        self.header_line = (
144
            f"UNITS LENGTH={self.units['length']} "
145
            f"ENERGY={self.units['energy']} "
146
            f"TIME={self.units['time']}"
147
        )
148

149 0
    def dump_to_file(self):
150

151 0
        self._initialize()
152

153
        # Loop over APR windows
154 0
        for windows in self.window_list:
155

156 0
            window, phase = parse_window(windows)
157

158
            # Check if file exist and write header line
159 0
            with open(os.path.join(self.path, windows, self.file_name), "w") as file:
160 0
                file.write(self.header_line + "\n")
161

162 0
            cv_index = 1
163 0
            cv_dict = {}
164 0
            cv_lines = []
165 0
            bias_lines = []
166

167 0
            self.group_index = 1
168 0
            self.group_atoms = {}
169

170
            # Parse each restraint in the list
171 0
            for restraint in self.restraint_list:
172
                # Skip restraint if the target or force constant is not defined.
173
                # Example: wall restraints only used during the attach phase.
174 0
                try:
175 0
                    target = restraint.phase[phase]["targets"][window]
176 0
                    force_constant = (
177
                        restraint.phase[phase]["force_constants"][window]
178
                        * self.k_factor
179
                    )
180 0
                except TypeError:
181 0
                    continue
182

183
                # Convert list to comma-separated string
184 0
                atom_index = self._get_atom_indices(restraint)
185 0
                atom_string = ",".join(map(str, atom_index))
186

187
                # Determine collective variable type
188 0
                colvar_type = "distance"
189 0
                if len(atom_index) == 3:
190 0
                    colvar_type = "angle"
191 0
                    target *= _PI_ / 180.0
192 0
                elif len(atom_index) == 4:
193 0
                    colvar_type = "torsion"
194 0
                    target *= _PI_ / 180.0
195

196
                # Determine bias type for this restraint
197 0
                bias_type = get_bias_potential_type(restraint, phase, window)
198

199
                # Append cv strings to lists
200
                # The code below prevents duplicate cv definition.
201
                # While not necessary, it makes the plumed file cleaner.
202 0
                if not get_key(cv_dict, atom_string):
203 0
                    cv_key = f"c{cv_index}"
204 0
                    cv_dict[cv_key] = atom_string
205

206 0
                    cv_lines.append(
207
                        f"{cv_key}: {colvar_type.upper()} ATOMS={atom_string} NOPBC\n"
208
                    )
209 0
                    bias_lines.append(
210
                        f"{bias_type.upper()} ARG={cv_key} AT={target:.4f} KAPPA={force_constant:.2f}\n"
211
                    )
212
                else:
213 0
                    cv_key = get_key(cv_dict, atom_string)[0]
214

215 0
                    bias_lines.append(
216
                        f"{bias_type.upper()} ARG={cv_key} AT={target:.4f} KAPPA={force_constant:.2f}\n"
217
                    )
218

219
                # Increment cv index
220 0
                cv_index += 1
221

222
            # Write collective variables to file
223 0
            self._write_colvar_to_file(windows, cv_lines, bias_lines)
224

225 0
    def _write_colvar_to_file(self, window, cv_list, bias_list):
226 0
        with open(os.path.join(self.path, window, self.file_name), "a") as file:
227 0
            if len(self.group_atoms) != 0:
228 0
                file.write("# Centroid groups\n")
229 0
                for key, value in self.group_atoms.items():
230 0
                    file.write(f"{key}: COM ATOMS={value}\n")
231

232 0
            file.write("# Collective variables\n")
233 0
            for line in cv_list:
234 0
                file.write(line)
235

236 0
            file.write("# Bias potentials\n")
237 0
            for line in bias_list:
238 0
                file.write(line)
239

240 0
    def _get_atom_indices(self, restraint):
241
        # Check atom index setting
242 0
        index_shift = 0
243 0
        if not restraint.amber_index:
244 0
            index_shift = 1
245 0
            logger.debug("Atom indices starts from 0 --> shifting indices by 1.")
246

247
        # Collect DAT atom indices
248 0
        atom_index = []
249

250 0
        if not restraint.group1:
251 0
            atom_index.append(restraint.index1[0] + index_shift)
252
        else:
253 0
            igr1 = ""
254 0
            for index in restraint.index1:
255 0
                igr1 += "{},".format(index + index_shift)
256

257 0
            if not get_key(self.group_atoms, igr1):
258 0
                self.group_atoms[f"g{self.group_index}"] = igr1
259 0
                self.group_index += 1
260

261 0
            atom_index.append(get_key(self.group_atoms, igr1)[0])
262

263 0
        if not restraint.group2:
264 0
            atom_index.append(restraint.index2[0] + index_shift)
265
        else:
266 0
            igr2 = ""
267 0
            for index in restraint.index2:
268 0
                igr2 += "{},".format(index + index_shift)
269

270 0
            if not get_key(self.group_atoms, igr2):
271 0
                self.group_atoms[f"g{self.group_index}"] = igr2
272 0
                self.group_index += 1
273

274 0
            atom_index.append(get_key(self.group_atoms, igr2)[0])
275

276 0
        if restraint.index3 and not restraint.group3:
277 0
            atom_index.append(restraint.index3[0] + index_shift)
278 0
        elif restraint.group3:
279 0
            igr3 = ""
280 0
            for index in restraint.index3:
281 0
                igr3 += "{},".format(index + index_shift)
282

283 0
            if not get_key(self.group_atoms, igr3):
284 0
                self.group_atoms[f"g{self.group_index}"] = igr3
285 0
                self.group_index += 1
286

287 0
            atom_index.append(get_key(self.group_atoms, igr3)[0])
288

289 0
        if restraint.index4 and not restraint.group4:
290 0
            atom_index.append(restraint.index4[0] + index_shift)
291 0
        elif restraint.group4:
292 0
            igr4 = ""
293 0
            for index in restraint.index4:
294 0
                igr4 += "{},".format(index + index_shift)
295

296 0
            if not get_key(self.group_atoms, igr4):
297 0
                self.group_atoms[f"g{self.group_index}"] = igr4
298 0
                self.group_index += 1
299

300 0
            atom_index.append(get_key(self.group_atoms, igr4)[0])
301

302 0
        return atom_index
303

304 0
    def add_dummy_atoms_to_file(self, structure):
305
        # Extract dummy atoms
306 0
        for windows in self.window_list:
307 0
            if isinstance(structure, str):
308 0
                structure = return_parmed_structure(structure)
309 0
            elif isinstance(structure, ParmedStructureClass):
310 0
                pass
311
            else:
312 0
                raise Exception(
313
                    "add_dummy_atoms_to_file does not support the type associated with structure: "
314
                    + type(structure)
315
                )
316

317 0
            dummy_atoms = extract_dummy_atoms(structure, serial=True)
318

319
            # Write dummy atom info to plumed file
320 0
            plumed_file = os.path.join(self.path, windows, self.file_name)
321 0
            if os.path.isfile(plumed_file):
322 0
                with open(plumed_file, "a") as file:
323 0
                    _write_dummy_to_file(file, dummy_atoms)
324
            else:
325 0
                raise Exception(f"ERROR: '{plumed_file}' file does not exists!")
326

327

328 0
def _check_plumed_units(units):
329
    """
330
    Checks the specified units and makes sure that its supported.
331
    """
332 0
    if units["energy"] not in ["kj/mol", "kcal/mol"]:
333 0
        raise Exception(
334
            f"Specified unit for energy ({units['energy']}) is not supported."
335
        )
336

337 0
    if units["length"] not in ["nm", "A"]:
338 0
        raise Exception(
339
            f"Specified unit for length ({units['length']}) is not supported."
340
        )
341

342 0
    if units["time"] not in ["ps", "fs", "ns"]:
343 0
        raise Exception(f"Specified unit for time ({units['time']}) is not supported.")
344

345

346 0
def _write_dummy_to_file(file, dummy_atoms, kpos=100.0):
347
    """
348
    Append to the "plumed.dat" file the dummy atoms colvar definition and position restraints
349

350
    Parameters
351
    ----------
352
    file : class '_io.TextIOWrapper'
353
        The file object handle to save the plumed file.
354
    dummy_atoms : dict
355
        Dictionary containing information about the dummy atoms.
356
    kpos : float
357
        Spring constant used to restrain dummy atoms (kcal/mol/A^2).
358

359
    Output
360
    ------
361
    dm1: POSITION ATOM = 170 NOPBC
362
    dm2: POSITION ATOM = 171 NOPBC
363
    dm3: POSITION ATOM = 172 NOPBC
364
    RESTRAINT...
365
    ARG = dm1.x, dm1.y, dm1.z, dm2.x, dm2.y, dm2.z, dm3.x, dm3.y, dm3.z
366
    AT = 19.68, 20.3, 26.9, 19.68, 20.3, 23.9, 19.68, 22.5, 21.7
367
    KAPPA = 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0
368
    ...
369
    RESTRAINT
370

371
    """
372

373 0
    file.write("# Dummy Atoms\n")
374 0
    file.write(f"dm1: POSITION ATOM={dummy_atoms['DM1']['idx']} NOPBC\n")
375 0
    file.write(f"dm2: POSITION ATOM={dummy_atoms['DM2']['idx']} NOPBC\n")
376 0
    file.write(f"dm3: POSITION ATOM={dummy_atoms['DM3']['idx']} NOPBC\n")
377

378 0
    arg = "dm1.x,dm1.y,dm1.z," "dm2.x,dm2.y,dm2.z," "dm3.x,dm3.y,dm3.z,"
379

380 0
    at = (
381
        f"{dummy_atoms['DM1']['pos'][0]:0.3f},"
382
        f"{dummy_atoms['DM1']['pos'][1]:0.3f},"
383
        f"{dummy_atoms['DM1']['pos'][2]:0.3f},"
384
    )
385 0
    at += (
386
        f"{dummy_atoms['DM2']['pos'][0]:0.3f},"
387
        f"{dummy_atoms['DM2']['pos'][1]:0.3f},"
388
        f"{dummy_atoms['DM2']['pos'][2]:0.3f},"
389
    )
390 0
    at += (
391
        f"{dummy_atoms['DM3']['pos'][0]:0.3f},"
392
        f"{dummy_atoms['DM3']['pos'][1]:0.3f},"
393
        f"{dummy_atoms['DM3']['pos'][2]:0.3f},"
394
    )
395

396 0
    kappa = (
397
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
398
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
399
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
400
    )
401

402 0
    file.write("RESTRAINT ...\n")
403 0
    file.write(f"ARG={arg}\n")
404 0
    file.write(f"AT={at}\n")
405 0
    file.write(f"KAPPA={kappa}\n")
406 0
    file.write("LABEL=dummy\n")
407 0
    file.write("... RESTRAINT\n")

Read our documentation on viewing source code .

Loading