readme.txt

#######################################################################################################################
##################################### This script was generated by MDPerTool v0.1 #####################################
#######################################################################################################################

import sys
import os
sys.path.append(r"C:\Users\law5_\anaconda3\envs\perto\Lib\site-packages\mdpertool")

from no_gui.get_positions_from_trajectory_file import *
from no_gui.energy_decomposition_from_trajectory import *
from no_gui.response_time_creator import *
from no_gui.Velocity_Changer import *
from no_gui.response_time_creator import *
from no_gui.apply_pdbfixer import fix_pdb
from no_gui.write_outputs import *
from no_gui.energy_component import *

from openmm import unit
from openmm import *
from openmm import app
from openmm.app import PME, NoCutoff, Ewald, CutoffPeriodic, CutoffNonPeriodic, HBonds, HAngles, AllBonds,DCDReporter, StateDataReporter
import openmm as mm
from openmm.unit import femtosecond, picosecond, nanometer, kelvin, angstrom, atmospheres, molar

from sys import stdout
import platform as pl
import time
import argparse
import multiprocessing as mp
from pathlib import Path, PureWindowsPath
from mdtraj import reporters
from mdtraj.reporters import XTCReporter
from src.logger import Logger
from datetime import datetime
from src.CustomReporters import Perturb_StateDataReporter

#######################################################################################################################
###################################### CLASSICAL MD PROCESS USING MDPerTool v0.1 ######################################
#######################################################################################################################

class Minimization_Reporter(openmm.MinimizationReporter):
    energies = []  # array to record progress
    logger_object = object

    def __init__(self, logger_object, *args):
        super().__init__(*args)
        self.energies = []  # array to record progress
        self.logger_object = logger_object

    def report(self, iteration, x, grad, args):
        # print current system energy to screen
        self.logger_object.info("Iteration: {} - Energy (kJ/mole): {}".format(iteration, args['system energy']))

        # save energy at each iteration to an array we can use later
        #self.energies.append(args['system energy'])

        # The report method must return a bool specifying if minimization should be stopped.
        # You can use this functionality for early termination.
        return False




# ---------------- LOGGER ---------------- #
current_date = datetime.now().strftime("%Y-%m-%d")
file_name = f"out_{current_date}.log"
log_path = os.path.join('C:/Users/law5_/Desktop/MDPerTool_GUI/Download', file_name)

log_obj = Logger(log_path)
"""
log_obj.critical("CRIC".format())
log_obj.error("ERR".format())
log_obj.warning("WARN".format())
log_obj.debug("debug".format())
log_obj.info("qwerty".format())
log_obj.info("asdfghjkl".format())
log_obj.info("zxcvbnm".format())
# closing file
log_obj.handlers.clear()
"""
# ---------------- LOGGER ---------------- #


simulation_last_time = 0
log_obj.info("pdb file fixing and preparing for simulation ...".format())
fixed_pdb_name = fix_pdb('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/5dk3_fixed_ph7.4.pdb', fixed_pdb_out_path='C:/Users/law5_/Desktop/MDPerTool_GUI/Download', logger_object=log_obj)

log_obj.info("Loading pdb to simulation engine ...".format())
pdb = app.PDBFile(fixed_pdb_name)

box = pdb.topology.getUnitCellDimensions()

log_obj.info("Modeller of pdb file is preparing ...".format())
modeller = mm.app.Modeller(pdb.topology, pdb.positions)
modeller.topology.setUnitCellDimensions(box)

log_obj.info("Forcefield parameters loading to the simulation system ...".format())
forcefield = app.ForceField('amber03.xml', 'tip3p.xml')

log_obj.critical("Adding missing hydrogens to the model ...".format())
modeller.addHydrogens(forcefield)

log_obj.info("Adding solvent (both water and ions) to the model to fill a rectangular box ...".format())
modeller.addSolvent(forcefield, model='tip3p', padding=10*angstrom, ionicStrength=0.1*molar)

log_obj.info("Constructing an OpenMM System ...".format())
system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,
                                      nonbondedCutoff=1.2*nanometer,
                                      constraints=None, rigidWater=True,
                                      ewaldErrorTolerance=0.005)

system.addForce(mm.MonteCarloBarostat(1 * atmospheres, 310.0*kelvin, 25))

nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
nonbonded.setUseSwitchingFunction(use=True)
nonbonded.setSwitchingDistance(1.0*nanometer)
nonbonded.setUseDispersionCorrection(True)

log_obj.info('Creating a %sIntegrator with %s %s . time step' %('Langevin', 2.0, 'femtosecond'))
integrator = mm.LangevinIntegrator(310.0*kelvin, 91.0/picosecond, 2.0*femtosecond)

if True == True:
    platform = mm.Platform.getPlatformByName('CUDA')
    properties = {'CudaPrecision': 'single'}
else:
    platform = mm.Platform.getPlatformByName('CUDA')
    properties = {'CUDAPrecision': 'single'}



simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)
simulation.context.computeVirtualSites()

log_obj.info('Minimizing for %s steps ...' % 500)
min_reporter = Minimization_Reporter(logger_object=log_obj)
simulation.minimizeEnergy(maxIterations=int(500), reporter=min_reporter)
log_obj.info("Minimization done, the energy is %s" % simulation.context.getState(getEnergy=True).getPotentialEnergy())
positions = simulation.context.getState(getPositions=True).getPositions()
log_obj.info("Minimized geometry is written to 'minimized.pdb'")
app.PDBFile.writeModel(modeller.topology, positions, open('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/minimized.pdb', 'w'), keepIds=True)

simulation.context.setVelocitiesToTemperature(310.0*kelvin)

log_obj.info('Equilibrating for %s steps ...' % 500)
simulation.step(int(500))

simulation.currentStep = 0
simulation.context.setTime(0)

log_obj.info('The trajectories will be saved in DCD file format.')
simulation.reporters.append(DCDReporter('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/output.dcd', 100))
log_obj.info("Saving DCD File for every 100 period")

log_obj.info("Saving XTC File for every 100 period")

log_obj.info("State Report will tell you.")
simulation.reporters.append(Perturb_StateDataReporter(stdout, 100, step=True, time=True,
                                                      potentialEnergy=True, kineticEnergy=True, totalEnergy=True,
                                                      temperature=True, progress=True, remainingTime=True, speed=True,
                                                      volume=True, density=True, totalSteps=1000,
                                                      report_type='Progress (%)'))

log_obj.info("Running Production...")
simulation.step(1000)
log_obj.info("Done!".format())

lastpositions = simulation.context.getState(getPositions=True).getPositions()

last_pdb = app.PDBFile.writeFile(modeller.topology, lastpositions, open('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/last.pdb', 'w'), keepIds=True)


state = simulation.context.getState(getPositions=True, getVelocities=True)

with open('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/system.xml', 'w') as f:
    system_xml = mm.XmlSerializer.serialize(system)
    f.write(system_xml)

with open('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/integrator.xml', 'w') as f:
    integrator_xml = mm.XmlSerializer.serialize(integrator)
    f.write(integrator_xml)

with open('C:/Users/law5_/Desktop/MDPerTool_GUI/Download/state.xml', 'w') as f:
    f.write(mm.XmlSerializer.serialize(state))

try:
    for reporter in simulation.reporters[:]:
        print("REPORTER: ", reporter)
        simulation.reporters.remove(reporter)
except Exception as E:
    print("EEEERRRORRRR: ", E)

#######################################################################################################################
#################################### PERTURBATION MD PROCESS USING MDPerTool v0.1 #####################################
#######################################################################################################################
# ## --> VARIABLES
simulation_last_step = 0
state_file_name = 'state.xml'
last_pdb = 'last.pdb'
dissipated_trajectory_name = 'energy_perturbation_trajectory'
undissipated_trajectory_name = 'without_energy_perturbation_trajectory'
reference_traj_file_for_pos = str
dissipation_traj_file_for_pos = str
OUTPUT_DIRECTORY = str
created_file_for_work = str
OUTPUT_FOLDER_NAME = str


OUTPUT_DIRECTORY = Path('C:/Users/law5_/Desktop/MDPerTool_GUI/Download')

last_pdb_file_path = os.path.join(OUTPUT_DIRECTORY, last_pdb)
modify_atoms = convert_res_to_atoms(last_pdb_file_path, ['GLN6'], 'CA')
state_file_path = os.path.join(OUTPUT_DIRECTORY, state_file_name)

for i in range(len([2])):
    name_of_changed_state_xml = change_velocity(state_file_path, [2][i], modify_atoms)
    new_dissipated_trajectory_name = dissipated_trajectory_name + str([2][i])

    if True == False and False == False:
        write_dcd_cond = True

    if i == 0:

        ################################################################################################################
        ################################ REFERENCE MD PROCESS USING MDPerTool v0.1 #####################################
        ################################################################################################################
        if 'CUDA' == 'OpenCL' and False == True:
            properties = {'OpenCLPrecision': 'double', 'OpenCLDeviceIndex': '1'}
            precision = 'double'

        if 'CUDA' == 'OpenCL' and False == False:
            properties = {'OpenCLPrecision': 'double'}
            precision = 'double'

        if 'CUDA' == 'CUDA' and False == True:
            properties = {'CudaPrecision': 'double', 'CudaDeviceIndex': '1'}
            precision = 'double'

        if 'CUDA' == 'CUDA' and False == False:
            properties = {'CudaPrecision': 'double'}
            precision = 'double'

        if 'CUDA' == 'CPU' and False == True:
            log_obj.info("The CPU platform always uses 'mixed' precision.".format())
            log_obj.info("Simulation process will use %s Thread(s)" % 2)
            properties = {'CpuThreads': '2'}
            precision = 'mixed'

        if 'CUDA' == 'Reference':
            log_obj.info("The Reference platform always uses 'double' precision.".format())
            properties = None
            precision= 'double'

        # we'll just take the topology from here...
        pdb = app.PDBFile(last_pdb_file_path)
        topology = pdb.topology

        log_obj.info("Forcefield parameters loading to the simulation system for NVE simulation ...".format())
        forcefield = app.ForceField('amber03.xml', 'tip3p.xml')

        log_obj.info("Constructing NVE System".format())
        reference_system = forcefield.createSystem(topology, nonbondedMethod=app.PME, nonbondedCutoff=1.2*nanometer,
                                                   constraints=None, rigidWater=True, ewaldErrorTolerance=1e-5)

        if True == True:
            log_obj.info("System will use Switching Distance".format())
            nonbonded = [f for f in reference_system.getForces() if isinstance(f, NonbondedForce)][0]
            nonbonded.setUseSwitchingFunction(use=True)
            nonbonded.setSwitchingDistance(1.0*nanometer)

        integrator = VerletIntegrator(1.0*femtosecond)
        integrator.setConstraintTolerance(1e-8)

        # let's specify our simulation platform again
        platform = mm.Platform.getPlatformByName('CUDA')

        # ok now let's do some simulation using this restraint
        if properties is None:
            ref_simulation = app.Simulation(topology, reference_system, integrator, platform)

        if properties is not None:
            ref_simulation = app.Simulation(topology, reference_system, integrator, platform, properties)

        ref_simulation.loadState(state_file_path)
        final_state = ref_simulation.context.getState(getVelocities=True, getPositions=True)
        positions = final_state.getPositions()
        velocities = final_state.getVelocities()
        ref_simulation.context.setPositions(positions)
        ref_simulation.context.computeVirtualSites()
        ref_simulation.context.setVelocities(velocities)


        ref_simulation.currentStep = simulation_last_step
        #ref_simulation.context.setTime(0)
        """
        if pl.system() == 'Windows':
            log_obj.info("Saving XTC File for every 1 period".format())
            XTC_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.xtc' % undissipated_trajectory_name)
            ref_simulation.reporters.append(XTCReporter(XTC_file_path, 1))
        """
        try:
            if True == False and False == False:
                log_obj.info("Saving DCD File for every 1 period".format())
                DCD_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.dcd' % undissipated_trajectory_name)
                ref_simulation.reporters.append(app.DCDReporter(DCD_file_path, 1))

            if True == True:
                log_obj.info("Saving DCD File for every 1 period".format())
                DCD_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.dcd' % undissipated_trajectory_name)
                ref_simulation.reporters.append(app.DCDReporter(DCD_file_path, 1))

            if False == True:
                log_obj.info("Saving XTC File for every 1 period".format())
                XTC_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.xtc' % undissipated_trajectory_name)
                ref_simulation.reporters.append(XTCReporter(XTC_file_path, 1))

        except Exception as e:
            log_obj.error(e, stack_info=True, exc_info=True)

        log_obj.info("State Report will tell you ...".format())
        ref_simulation.reporters.append(Perturb_StateDataReporter(stdout, 1, step=True, time=True, potentialEnergy=True,
                                  kineticEnergy=True, totalEnergy=True, temperature=True, progress=True, volume=True,
                                  density=True, remainingTime=True, speed=True, totalSteps=1000,
                                  report_type='Reference MD Progress (%)'))

        ref_simulation.step(1000)

        simulation_last_step = 0

    try:
        for reporter in ref_simulation.reporters[:]:
            print("REPORTER: ", reporter)
            ref_simulation.reporters.remove(reporter)
    except Exception as E:
        print("EEEERRRORRRR: ", E)

    ################################################################################################################
    ############################### DISSIPATION MD PROCESS USING MDPerTool v0.1 ####################################
    ################################################################################################################

    if 'CUDA' == 'OpenCL' and False == True:
        properties = {'OpenCLPrecision': 'double', 'OpenCLDeviceIndex': '1'}
        precision = 'double'

    if 'CUDA' == 'OpenCL' and False == False:
        properties = {'OpenCLPrecision': 'double'}
        precision = 'double'

    if 'CUDA' == 'CUDA' and False == True:
        properties = {'CudaPrecision': 'double', 'CudaDeviceIndex': '1'}
        precision = 'double'

    if 'CUDA' == 'CUDA' and False == False:
        properties = {'CudaPrecision': 'double'}
        precision = 'double'

    if 'CUDA' == 'CPU' and False == True:
        log_obj.info("The CPU platform always uses 'mixed' precision.".format())
        log_obj.info("Simulation process will use %s Thread(s)" % 2)
        properties = {'CpuThreads': '2'}
        precision = 'mixed'

    if 'CUDA' == 'Reference':
        log_obj.info("The Reference platform always uses 'double' precision.")
        properties = None
        precision= 'double'

    log_obj.info("System will use %s Platform with %s Precision" % ('CUDA', precision))

    # we'll just take the topology from here...
    pdb = app.PDBFile(last_pdb_file_path)
    topology = pdb.topology

    log_obj.info("Forcefield parameters loading to the simulation system ...".format())
    forcefield = app.ForceField('amber03.xml', 'tip3p.xml')

    log_obj.info("Constructing an OpenMM System".format())
    perturbed_system = forcefield.createSystem(topology, nonbondedMethod=app.PME, nonbondedCutoff=1.2*nanometer,
                                               constraints=None, rigidWater=True, ewaldErrorTolerance=1e-5)

    if True == True:
        log_obj.info("System will use Switching Distance".format())
        nonbonded = [f for f in perturbed_system.getForces() if isinstance(f, NonbondedForce)][0]
        nonbonded.setUseSwitchingFunction(use=True)
        nonbonded.setSwitchingDistance(1.0*nanometer)

    integrator = VerletIntegrator(1.0*femtosecond)
    integrator.setConstraintTolerance(1e-8)

    # let's specify our simulation platform again
    platform = mm.Platform.getPlatformByName('CUDA')

    # ok now let's do some simulation using this restraint
    if properties is None:
        dis_simulation = app.Simulation(topology, perturbed_system, integrator, platform)

    if properties is not None:
        dis_simulation = app.Simulation(topology, perturbed_system, integrator, platform, properties)

    dis_simulation.loadState(name_of_changed_state_xml)
    final_state = dis_simulation.context.getState(getVelocities=True, getPositions=True)
    positions = final_state.getPositions()
    velocities = final_state.getVelocities()
    dis_simulation.context.setPositions(positions)
    dis_simulation.context.computeVirtualSites()
    dis_simulation.context.setVelocities(velocities)
    log_obj.info("Will the perturbed system use periodic boundary conditions? --> %s" % perturbed_system.usesPeriodicBoundaryConditions())

    dis_simulation.currentStep = 0
    dis_simulation.context.setTime(0)

    """
    if pl.system() == 'Windows':
        log_obj.info("Saving XTC File for every 1 period")
        XTC_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.xtc' % new_dissipated_trajectory_name)
        dis_simulation.reporters.append(XTCReporter(XTC_file_path, 1))
    """

    try:
        if True == False and False == False:
            log_obj.info("Saving DCD File for every 1 period")
            DCD_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.dcd' % new_dissipated_trajectory_name)
            dis_simulation.reporters.append(app.DCDReporter(DCD_file_path, 1))

        if True == True:
            log_obj.info("Saving DCD File for every 1 period")
            DCD_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.dcd' % new_dissipated_trajectory_name)
            dis_simulation.reporters.append(app.DCDReporter(DCD_file_path, 1))

        if False == True:
            log_obj.info("Saving XTC File for every 1 period")
            XTC_file_path = os.path.join(OUTPUT_DIRECTORY, '%s.xtc' % new_dissipated_trajectory_name)
            dis_simulation.reporters.append(XTCReporter(XTC_file_path, 1))

    except Exception as e:
        log_obj.error(e, stack_info=True, exc_info=True)

    log_obj.info("State Report will tell you ...")
    dis_simulation.reporters.append(Perturb_StateDataReporter(stdout, 1, step=True, time=True, potentialEnergy=True,
                              kineticEnergy=True, totalEnergy=True, temperature=True, progress=True, volume=True,
                              density=True, remainingTime=True, speed=True, totalSteps=1000,
                              report_type='Dissipation MD Progress (%)'))

    dis_simulation.step(1000)

    simulation_last_step = 0

    try:
        for reporter in dis_simulation.reporters[:]:
            print("REPORTER: ", reporter)
            dis_simulation.reporters.remove(reporter)
    except Exception as E:
        print("EEEERRRORRRR: ", E)

    del dis_simulation, ref_simulation
    ####################################################################################################################
    ############################### PER RESIDUE ENERGY CALCULATION USING MDPerTool v0.1 ################################
    ####################################################################################################################
    """
    if pl.system() == 'Windows':
        log_obj.info("Decompose started using XTC File ...".format())
        reference_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, undissipated_trajectory_name + '.xtc')
        dissipation_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, new_dissipated_trajectory_name + '.xtc')
    """
    try:
        if True == False and False == False:
            log_obj.info("Decompose started using DCD File ...".format())
            reference_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, undissipated_trajectory_name + '.dcd').replace('\\','/')
            dissipation_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, new_dissipated_trajectory_name + '.dcd').replace('\\','/')

        if True == True:
            log_obj.info("Decompose started using DCD File ....".format())
            reference_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, undissipated_trajectory_name + '.dcd').replace('\\','/')
            dissipation_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, new_dissipated_trajectory_name + '.dcd').replace('\\','/')
            print(reference_traj_file_for_pos)
            print(dissipation_traj_file_for_pos)
            print(last_pdb_file_path)

        if False == True:
            log_obj.info("Decompose started using XTC File ....".format())
            reference_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, undissipated_trajectory_name + '.xtc').replace('\\','/')
            dissipation_traj_file_for_pos = os.path.join(OUTPUT_DIRECTORY, new_dissipated_trajectory_name + '.xtc').replace('\\','/')

    except Exception as e:
        print("ERROR: ", e)
        log_obj.error(e, stack_info=True, exc_info=True)

    try:

        #if i == 0:
        #    position_list, unwrap_pdb = get_openmm_pos_from_traj_with_mdtraj(top=last_pdb_file_path,
        #                                                                     ref_traj=reference_traj_file_for_pos,
        #                                                                     modif_traj=dissipation_traj_file_for_pos,
        #                                                                     logger_object=log_obj)


        #if i != 0:
        #    position_list, unwrap_pdb = get_openmm_pos_from_traj_with_mdtraj(top=last_pdb_file_path,
        #                                                             ref_traj=None,
        #                                                             modif_traj=dissipation_traj_file_for_pos,
        #                                                             logger_object=log_obj)

         # --> RESIDUE BASED DECOMPOSITION
        if i == 0:

            #residue_based_decomposition(topol=unwrap_pdb, trj_pos_list=position_list, start_res=0, stop_res=250,
            #                            output_directory=OUTPUT_DIRECTORY, que=None, platform_name='CUDA',
            #                            ref_energy_name='reference_energy_file.csv', device_id_active=False,
            #                            num_of_threads=2,
            #                            modif_energy_name='modified_energy_file_%s.csv' % int([2][i]),
            #                            origin_last_pdb=last_pdb_file_path, ff='amber03.xml', logger_object=log_obj)

            process_energy_data(topology_file=last_pdb_file_path, protein_ff='amber03.xml', water_ff='tip3p.xml',
                                modif_trajectory=dissipation_traj_file_for_pos, device_id_active=False,
                                ref_trajectory=reference_traj_file_for_pos, platform_type='CUDA',
                                nonbonded_cutoff=1.2*nanometer, ref_energy_name = 'reference_energy_file.csv',
                                modif_energy_name='modified_energy_file_%s.csv' % int([2][i]), 
                                output_directory=OUTPUT_DIRECTORY, num_of_threads=2, que=None, 
                                logger_object=log_obj)
	



        if i != 0:
            # --> RESIDUE BASED DECOMPOSITION
            """
            residue_based_decomposition(topol=unwrap_pdb, trj_pos_list=position_list, start_res=0, stop_res=250,
                                        output_directory=OUTPUT_DIRECTORY, que=None, platform_name='CUDA',
                                        ref_energy_name=None, device_id_active=False, num_of_threads=2,
                                        modif_energy_name='modified_energy_file_%s.csv' % int([2][i]),
                                        origin_last_pdb=last_pdb_file_path,
                                        ff='amber03.xml', logger_object=log_obj)

            """

            process_energy_data(topology_file=last_pdb_file_path, protein_ff='amber03.xml', water_ff='tip3p.xml',
                                modif_trajectory=dissipation_traj_file_for_pos, device_id_active=False,
                                ref_trajectory=reference_traj_file_for_pos, platform_type='CUDA',
                                nonbonded_cutoff=1.2*nanometer, ref_energy_name = 'reference_energy_file.csv',
                                modif_energy_name='modified_energy_file_%s.csv' % int([2][i]), 
                                output_directory=OUTPUT_DIRECTORY, num_of_threads=2, que=None, 
                                logger_object=log_obj)




        # --> RESPONSE TIME CSV EXPORTER
        getResidueResponseTimes(os.path.join(OUTPUT_DIRECTORY, 'reference_energy_file.csv'),
                                os.path.join(OUTPUT_DIRECTORY, 'modified_energy_file_%s.csv' % int([2][i])),
                                outputName=os.path.join(OUTPUT_DIRECTORY, 'responseTimes_%s.csv' % int([2][i])))

    except Exception as e:
        print("ERROR: ", e)
        log_obj.error(e, stack_info=True, exc_info=True)

    try:
        del XTC_file_path, ref_simulation, dis_simulation, simulation
    except NameError:
        pass