Python Molecule.atom_select方法代码示例

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import atom_select [as 别名]
    def __init__(self,options,tgt_opts,forcefield):
        super(BindingEnergy,self).__init__(options,tgt_opts,forcefield)
        self.set_option(None, None, 'inter_txt', os.path.join(self.tgtdir,tgt_opts['inter_txt']))
        self.global_opts, self.sys_opts, self.inter_opts = parse_interactions(self.inter_txt)
        # If the global option doesn't exist in the system / interaction, then it is copied over.
        for opt in self.global_opts:
            for sys in self.sys_opts:
                if opt not in self.sys_opts[sys]:
                    self.sys_opts[sys][opt] = self.global_opts[opt]
            for inter in self.inter_opts:
                if opt not in self.inter_opts[inter]:
                    self.inter_opts[inter][opt] = self.global_opts[opt]
        for inter in self.inter_opts:
            if 'energy_unit' in self.inter_opts[inter] and self.inter_opts[inter]['energy_unit'].lower() not in ['kilocalorie_per_mole', 'kilocalories_per_mole']:
                logger.error('Usage of physical units is has been removed, please provide all binding energies in kcal/mole\n')
                raise RuntimeError
            self.inter_opts[inter]['reference_physical'] = self.inter_opts[inter]['energy']

        if tgt_opts['energy_denom'] == 0.0:
            self.set_option(None, None, 'energy_denom', val=np.std(np.array([val['reference_physical'] for val in self.inter_opts.values()])))
        else:
            self.set_option(None, None, 'energy_denom', val=tgt_opts['energy_denom'])

        self.set_option(None, None, 'rmsd_denom', val=tgt_opts['rmsd_denom'])

        self.set_option(tgt_opts,'cauchy')
        self.set_option(tgt_opts,'attenuate')
        ## LPW 2018-02-11: This is set to True if the target calculates
        ## a single-point property over several existing snapshots.
        self.loop_over_snapshots = False

        logger.info("The energy denominator is: %s\n" % str(self.energy_denom)) 
        logger.info("The RMSD denominator is: %s\n" % str(self.rmsd_denom))

        if self.cauchy:
            logger.info("Each contribution to the interaction energy objective function will be scaled by 1.0 / ( denom**2 + reference**2 )\n")
            if self.attenuate:
                logger.error('attenuate and cauchy are mutually exclusive\n')
                raise RuntimeError
        elif self.attenuate:
            logger.info("Repulsive interactions beyond energy_denom will be scaled by 1.0 / ( denom**2 + (reference-denom)**2 )\n")
        ## Build keyword dictionaries to pass to engine.
        engine_args = OrderedDict(list(self.OptionDict.items()) + list(options.items()))
        del engine_args['name']
        ## Create engine objects.
        self.engines = OrderedDict()
        for sysname,sysopt in self.sys_opts.items():
            M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']))
            if 'select' in sysopt:
                atomselect = np.array(uncommadash(sysopt['select']))
                M = M.atom_select(atomselect)
            if self.FF.rigid_water: M.rigid_water()
            self.engines[sysname] = self.engine_(target=self, mol=M, name=sysname, tinker_key=os.path.join(sysopt['keyfile']), **engine_args)

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import atom_select [as 别名]
 def prepare_temp_directory(self, options, tgt_opts):
     abstempdir = os.path.join(self.root,self.tempdir)
     if self.FF.rigid_water:
         self.optprog = "optrigid"
         #LinkFile(os.path.join(self.root,self.tgtdir,"rigid.key"),os.path.join(abstempdir,"rigid.key"))
     else:
         self.optprog = "optimize"
     # Link the necessary programs into the temporary directory
     LinkFile(os.path.join(options['tinkerpath'],"analyze"),os.path.join(abstempdir,"analyze"))
     LinkFile(os.path.join(options['tinkerpath'],self.optprog),os.path.join(abstempdir,self.optprog))
     LinkFile(os.path.join(options['tinkerpath'],"superpose"),os.path.join(abstempdir,"superpose"))
     # Link the run parameter file
     # The master file might be unneeded??
     for sysname,sysopt in self.sys_opts.items():
         if self.FF.rigid_water:
             # Make every water molecule rigid
             # WARNING: Hard coded values here!
             M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']),ftype="tinker")
             for a in range(0, len(M.xyzs[0]), 3):
                 flex = M.xyzs[0]
                 wat = flex[a:a+3]
                 com = wat.mean(0)
                 wat -= com
                 o  = wat[0]
                 h1 = wat[1]
                 h2 = wat[2]
                 r1 = h1 - o
                 r2 = h2 - o
                 r1 /= Np.linalg.norm(r1)
                 r2 /= Np.linalg.norm(r2)
                 # Obtain unit vectors.
                 ex = r1 + r2
                 ey = r1 - r2
                 ex /= Np.linalg.norm(ex)
                 ey /= Np.linalg.norm(ey)
                 Bond = 0.9572
                 Ang = Np.pi * 104.52 / 2 / 180
                 cosx = Np.cos(Ang)
                 cosy = Np.sin(Ang)
                 h1 = o + Bond*ex*cosx + Bond*ey*cosy
                 h2 = o + Bond*ex*cosx - Bond*ey*cosy
                 rig = Np.array([o, h1, h2]) + com
                 M.xyzs[0][a:a+3] = rig
             M.write(os.path.join(abstempdir,sysopt['geometry']),ftype="tinker")
         else:
             M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']),ftype="tinker")
             if 'select' in sysopt:
                 atomselect = Np.array(uncommadash(sysopt['select']))
                 #atomselect = eval("Np.arange(M.na)"+sysopt['select'])
                 M = M.atom_select(atomselect)
             M.write(os.path.join(abstempdir,sysname+".xyz"),ftype="tinker")
             write_key_with_prm(os.path.join(self.root,self.tgtdir,sysopt['keyfile']),os.path.join(abstempdir,sysname+".key"),ffobj=self.FF)

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import atom_select [as 别名]

#.........这里部分代码省略.........
            # Two equivalent ways to get the RMSD, here for reference.
            #----
            # M1 = Molecule("%s.xyz" % self.name, ftype="tinker")
            # M2 = Molecule("%s.xyz_2" % self.name, ftype="tinker")
            # M1 += M2
            # rmsd = M1.ref_rmsd(0)[1]
            #----
            # oo = self.calltinker("superpose %s.xyz %s.xyz_2 1 y u n 0" % (self.name, self.name))
            # for line in oo:
            #     if "Root Mean Square Distance" in line:
            #         rmsd = float(line.split()[-1])
            #----
            os.system("rm %s.xyz_2" % self.name)
        else:
            o = self.calltinker("analyze %s.xyz E" % self.name)
        # Read the TINKER output. 
        E = None
        for line in o:
            if "Total Potential Energy" in line:
                E = float(line.split()[-2].replace('D','e'))
        if E == None:
            logger.error("Total potential energy wasn't encountered when calling analyze!\n")
            raise RuntimeError
        if optimize and abs(E-E_) > 0.1:
            warn_press_key("Energy from optimize and analyze aren't the same (%.3f vs. %.3f)" % (E, E_))
        return E, rmsd

    def interaction_energy(self, fraga, fragb):
        
        """ Calculate the interaction energy for two fragments. """

        logger.error('Interaction energy is not yet implemented in AMBER interface')
        raise NotImplementedError
        self.A = TINKER(name="A", mol=self.mol.atom_select(fraga), tinker_key="%s.key" % self.name, tinkerpath=self.tinkerpath)
        self.B = TINKER(name="B", mol=self.mol.atom_select(fragb), tinker_key="%s.key" % self.name, tinkerpath=self.tinkerpath)

        # Interaction energy needs to be in kcal/mol.
        return (self.energy() - self.A.energy() - self.B.energy()) / 4.184

    def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, nequil=0, nsave=1000, minimize=True, anisotropic=False, threads=1, verbose=False, **kwargs):
        
        """
        Method for running a molecular dynamics simulation.  

        Required arguments:
        nsteps      = (int)   Number of total time steps
        timestep    = (float) Time step in FEMTOSECONDS
        temperature = (float) Temperature control (Kelvin)
        pressure    = (float) Pressure control (atmospheres)
        nequil      = (int)   Number of additional time steps at the beginning for equilibration
        nsave       = (int)   Step interval for saving and printing data
        minimize    = (bool)  Perform an energy minimization prior to dynamics
        threads     = (int)   Specify how many OpenMP threads to use

        Returns simulation data:
        Rhos        = (array)     Density in kilogram m^-3
        Potentials  = (array)     Potential energies
        Kinetics    = (array)     Kinetic energies
        Volumes     = (array)     Box volumes
        Dips        = (3xN array) Dipole moments
        EComps      = (dict)      Energy components
        """

        logger.error('Molecular dynamics not yet implemented in AMBER interface')
        raise NotImplementedError

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import atom_select [as 别名]

#.........这里部分代码省略.........
            o = self.calltinker("analyze %s.xyz_2 E" % self.name)
            #----
            # Two equivalent ways to get the RMSD, here for reference.
            #----
            # M1 = Molecule("%s.xyz" % self.name, ftype="tinker")
            # M2 = Molecule("%s.xyz_2" % self.name, ftype="tinker")
            # M1 += M2
            # rmsd = M1.ref_rmsd(0)[1]
            #----
            # oo = self.calltinker("superpose %s.xyz %s.xyz_2 1 y u n 0" % (self.name, self.name))
            # for line in oo:
            #     if "Root Mean Square Distance" in line:
            #         rmsd = float(line.split()[-1])
            #----
            os.system("rm %s.xyz_2" % self.name)
        else:
            o = self.calltinker("analyze %s.xyz E" % self.name)
        # Read the TINKER output. 
        E = None
        for line in o:
            if "Total Potential Energy" in line:
                E = float(line.split()[-2].replace('D','e'))
        if E == None:
            logger.error("Total potential energy wasn't encountered when calling analyze!\n")
            raise RuntimeError
        if optimize and abs(E-E_) > 0.1:
            warn_press_key("Energy from optimize and analyze aren't the same (%.3f vs. %.3f)" % (E, E_))
        return E, rmsd

    def interaction_energy(self, fraga, fragb):
        
        """ Calculate the interaction energy for two fragments. """

        self.A = TINKER(name="A", mol=self.mol.atom_select(fraga), tinker_key="%s.key" % self.name, tinkerpath=self.tinkerpath)
        self.B = TINKER(name="B", mol=self.mol.atom_select(fragb), tinker_key="%s.key" % self.name, tinkerpath=self.tinkerpath)

        # Interaction energy needs to be in kcal/mol.
        return (self.energy() - self.A.energy() - self.B.energy()) / 4.184

    def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, nequil=0, nsave=1000, minimize=True, anisotropic=False, threads=1, verbose=False, **kwargs):
        
        """
        Method for running a molecular dynamics simulation.  

        Required arguments:
        nsteps      = (int)   Number of total time steps
        timestep    = (float) Time step in FEMTOSECONDS
        temperature = (float) Temperature control (Kelvin)
        pressure    = (float) Pressure control (atmospheres)
        nequil      = (int)   Number of additional time steps at the beginning for equilibration
        nsave       = (int)   Step interval for saving and printing data
        minimize    = (bool)  Perform an energy minimization prior to dynamics
        threads     = (int)   Specify how many OpenMP threads to use

        Returns simulation data:
        Rhos        = (array)     Density in kilogram m^-3
        Potentials  = (array)     Potential energies
        Kinetics    = (array)     Kinetic energies
        Volumes     = (array)     Box volumes
        Dips        = (3xN array) Dipole moments
        EComps      = (dict)      Energy components
        """

        md_defs = OrderedDict()
        md_opts = OrderedDict()
        # Print out averages only at the end.