Commit b4871d69 authored by milafternoon's avatar milafternoon
Browse files

fixing stuff to enable RMSE derivatives

parent 5c5c5cfd
Loading
Loading
Loading
Loading
+5 −1
Original line number Diff line number Diff line
@@ -207,7 +207,11 @@ class Mod:
        """
        for t in self.top:
            self.topname = path
            t.save_mod(self.topname + name + '-' + t.top)
            final_name = self.topname + name + '-' + t.top
            if not os.path.exists(final_name):
                t.save_mod(final_name)
            else:
                print(f'file {final_name} already exists, leaving as it is')
    
    def goto_mydir(self):
        """
+33 −29
Original line number Diff line number Diff line
@@ -23,16 +23,20 @@ class ModNbfix(Mod):
        self.nbe = 'm' if 'm' in changes else ''
        if not (self.nbs or self.nbe):
            raise ValueError("no mode selected appropriate for an ModParam object")
        for t in self.top:
        for t, tnames, tres, tresnames in zip(self.top, self.names, self.res, self.resnames):
            self.LJ_sigmas = []
            self.LJ_epsilons = []
            self.NB_sigmas = None
            self.NB_epsilons = None
            self.types_are_consistent(t)
            self.types_are_consistent(t, tnames, tres, tresnames)
            if 'nonbond_params' not in t.section_headers:
                self.add_nbfix_section(t)
            self.types = (self.get_type(t, self.names[0][0], self.res[0][0], self.resnames[0][0]),
                          self.get_type(t, self.names[1][0], self.res[1][0], self.resnames[1][0]))
            try:
                self.types = (self.get_type(t, tnames[0][0], tres[0][0], tresnames[0][0]),
                              self.get_type(t, tnames[1][0], tres[1][0], tresnames[1][0]))
            except IndexError:
                pass
            else:
                self.check_lj_params(t)
                if self.types[0] == self.types[1]:
                    self.prefixes = ['Y', 'Y']
@@ -40,11 +44,11 @@ class ModNbfix(Mod):
                    self.prefixes = ['Q', 'Y']
                # we mod the first type
                self.clone_type(t, atomtype=self.types[0], prefix=self.prefixes[0])
            self.prefix_type(t, 0)
                self.prefix_type(t, tnames, tres, tresnames, 0)
                # and then second, if different than first
                if self.types[0] != self.types[1]:
                    self.clone_type(t, atomtype=self.types[1], prefix=self.prefixes[1])
                self.prefix_type(t, 1)
                    self.prefix_type(t, tnames, tres, tresnames, 1)
                self.sort_dihedrals(t)
                self.mod_nb_sigma_eps(t)
    
@@ -100,16 +104,16 @@ class ModNbfix(Mod):
        returnlist.append(fstring.format(*mod_list))
        return returnlist
    
    def prefix_type(self, top, prefix_number):
    def prefix_type(self, top, tnames, tres, tresnames, prefix_number):
        """
        looks for a line in topology section 'atoms'
        and modifies it as needed
        :return: None
        """
        for x in range(len(self.names[prefix_number])):
            name = self.names[prefix_number][x]
            resname = self.resnames[prefix_number][x]
            resnum = self.res[prefix_number][x]
        for x in range(len(tnames[prefix_number])):
            name = tnames[prefix_number][x]
            resname = tresnames[prefix_number][x]
            resnum = tres[prefix_number][x]
            for section in top.list_sections('atoms'):
                for line in range(len(top.sections[section])):
                    lspl = top.sections[section][line].split()
@@ -141,18 +145,18 @@ class ModNbfix(Mod):
                raise RuntimeError("You need to provide mass for the atom to be changed")
            return fstring.format(*cont_list[:8], prefix + cont_list[1], cont_list[6], cont_list[7])
    
    def types_are_consistent(self, top):
    def types_are_consistent(self, top, tnames, tres, tresnames):
        """
        if many atoms are passed in a single input line, need to make sure
        they are all of the same type; raises ValueError if this is not the case
        :return: None
        """
        consistent1 = all([self.get_type(top, self.names[0][0], self.res[0][0], self.resnames[0][0])
                          == self.get_type(top, self.names[0][n], self.res[0][n], self.resnames[0][n])
                          for n in range(len(self.names[0]))])
        consistent2 = all([self.get_type(top, self.names[1][0], self.res[1][0], self.resnames[1][0])
                           == self.get_type(top, self.names[1][n], self.res[1][n], self.resnames[1][n])
                           for n in range(len(self.names[1]))])
        consistent1 = all([self.get_type(top, tnames[0][0], tres[0][0], tresnames[0][0])
                          == self.get_type(top, tnames[0][n], tres[0][n], tresnames[0][n])
                          for n in range(len(tnames[0]))])
        consistent2 = all([self.get_type(top, tnames[1][0], tres[1][0], tresnames[1][0])
                           == self.get_type(top, tnames[1][n], tres[1][n], tresnames[1][n])
                           for n in range(len(tnames[1]))])
        if not consistent1 or not consistent2:
            raise ValueError("atoms within a single nbfix modification need to have consistent types")

+20 −15
Original line number Diff line number Diff line
@@ -100,6 +100,9 @@ def run_gromacs(args):
                            '-dhdl {fname}/{tn}-dhdl >> rerun.log 2>&1'.format(gmx=gmx, fname=folder_name,
                                                                               traj=traj_path, tn=traj_alias,
                                                                               tpr=tpr), shell=True)
        else:
            print("in mod {fname} skipping trajectory {tn}, already calculated".format(fname=folder_name,
                                                                                       tn=traj_alias))
        
        
def get_pot_code(err_msg):
@@ -300,33 +303,35 @@ def parse_input_line(line, top, td):
        delim = ';'
        for mod in mods:
            if c == 0:
                t1 = line_list[2]
                t2 = line_list[3]
                ll1 = top.find_type(t1)
                ll2 = top.find_type(t2)
                n = [tuple(x.split('-')[2] for x in ll1), tuple(x.split('-')[2] for x in ll2)]
                rn = [tuple(int(x.split('-')[1]) for x in ll1), tuple(int(x.split('-')[1]) for x in ll2)]
                r = [tuple(x.split('-')[0] for x in ll1), tuple(x.split('-')[0] for x in ll2)]
                t1 = line_list[1]
                t2 = line_list[2]
                n, rn, r = [], [], []
                for t in top:
                    ll1 = t.find_type(t1)
                    ll2 = t.find_type(t2)
                    n.append([tuple(x.split('-')[2] for x in ll1), tuple(x.split('-')[2] for x in ll2)])
                    rn.append([tuple(int(x.split('-')[1]) for x in ll1), tuple(int(x.split('-')[1]) for x in ll2)])
                    r.append([tuple(x.split('-')[0] for x in ll1), tuple(x.split('-')[0] for x in ll2)])
            elif c == 1:
                assert line.count(delim) == 1
                assert line.count(delim) == 1 and len(top) == 1  # TODO not clear how to implement for more tops, think
                line_first, line_second = line.split(delim)
                line_first, line_second = line_first.split(), line_second.split()
                ll1 = line_first[2:]
                ll2 = line_second[:]
                num = int(line_list[1].split('-')[1])
                res = line_list[1].split('-')[0]
                n = [tuple(ll1), tuple(ll2)]
                rn = [tuple(num for _ in ll1), tuple(num for _ in ll2)]
                r = [tuple(res for _ in ll1), tuple(res for _ in ll2)]
                n = [[tuple(ll1), tuple(ll2)]]
                rn = [[tuple(num for _ in ll1), tuple(num for _ in ll2)]]
                r = [[tuple(res for _ in ll1), tuple(res for _ in ll2)]]
            elif c == 2:
                assert line.count(delim) == 1
                assert line.count(delim) == 1 and len(top) == 1
                line_first, line_second = line.split(delim)
                line_first, line_second = line_first.split(), line_second.split()
                ll1 = line_first[1:]
                ll2 = line_second[:]
                n = [tuple(x.split('-')[2] for x in ll1), tuple(x.split('-')[2] for x in ll2)]
                rn = [tuple(int(x.split('-')[1]) for x in ll1), tuple(x.split('-')[1] for x in ll2)]
                r = [tuple(x.split('-')[0] for x in ll1), tuple(x.split('-')[0] for x in ll2)]
                n = [[tuple(x.split('-')[2] for x in ll1), tuple(x.split('-')[2] for x in ll2)]]
                rn = [[tuple(int(x.split('-')[1]) for x in ll1), tuple(x.split('-')[1] for x in ll2)]]
                r = [[tuple(x.split('-')[0] for x in ll1), tuple(x.split('-')[0] for x in ll2)]]
            else:
                raise ValueError("input line {} cannot be processed".format(' '.join(line_list)))
            new_mods.append(ModNbfix(top=deepcopy(top), master=td, resnr=rn, names=n, resnames=r, changes=mod))
+8 −3
Original line number Diff line number Diff line
@@ -158,9 +158,14 @@ class Quantity:
                    self.discrete_free_energy_derivative = np.array(self.multi_discrete_free_derivative())

    def calc_discrete_derivs(self):
        mean_der = [np.sum(np.concatenate([traj.derivatives * (traj.rc == state) for traj in self.trajs])) /
                    np.sum(np.concatenate([traj.rc == state for traj in self.trajs])) for state in self.discrete_states]
        mean_der = [np.sum(np.concatenate([traj.derivatives * traj.weights * (traj.rc == state)
                                           for traj in self.trajs])) /
                    np.sum(np.concatenate([traj.weights * (traj.rc == state)
                                           for traj in self.trajs])) for state in self.discrete_states]
        if len(self.discrete_states) > 1:
            rel_derivatives = [mean_der[0] - der for der in mean_der[1:]]
        else:
            rel_derivatives = mean_der[0]
        return np.array([[deriv, 0] for deriv in rel_derivatives]).ravel()  # TODO implement stderr

    def calc_ens_avg_profile(self, attr_name, second_attr=None):
+1 −1
Original line number Diff line number Diff line
@@ -23,7 +23,7 @@ class ThermoDiff:
    """
    def __init__(self, options):
        self.wdir = os.getcwd()
        self.topol = [self.abs_path(x) for x in options.topol.split()]
        self.topol = [self.abs_path(x) for x in options.topol.split()]  # TODO should crash if topols have identical names
        self.struct = [self.abs_path(x) for x in options.struct.split()]
        self.list_mods = self.abs_path(options.list_mods)
        # need a bit of workaround in case trajs_list file is in a different directory than the actual trajectories
Loading