diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 0603e70..fec5f40 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -19,10 +19,10 @@ jobs: steps: - uses: actions/checkout@v4 - - name: Set up Python 3.10 + - name: Set up Python 3.13 uses: actions/setup-python@v3 with: - python-version: "3.10" + python-version: "3.13" - name: Install pip dependencies run: | python -m pip install --upgrade pip diff --git a/README.md b/README.md index 675f4a1..77d5d6e 100644 --- a/README.md +++ b/README.md @@ -26,14 +26,18 @@ The examples described in the paper can be found in the `examples` folder. ## Installation Instructions -1. Make new conda environment for calvados +1. Create a new conda environment for calvados ``` -conda create -n calvados python=3.10 +conda create -n calvados python=3.13 conda activate calvados ``` -(2. Only needed when planning to use GPUs: Install openmm via conda-force with cudatoolkit. This step can be skipped if running on CPU only.) +2. Install openmm via conda-force with cudatoolkit, if you are planning to run on GPUs with CUDA ``` -conda install -c conda-forge openmm=8.2.0 cudatoolkit=11.8 +conda install -c conda-forge openmm=8.2.0 cudatoolkit=11.8 mdanalysis=2.9 mdtraj=1.11 +``` +or without cudatoolkit if you are running on CPU or OpenCL +``` +conda install -c conda-forge openmm=8.2.0 mdanalysis=2.9 mdtraj=1.11 ``` 3. Clone package and install CALVADOS and its dependencies using pip ``` diff --git a/calvados/analysis.py b/calvados/analysis.py index 8a7df19..19af9ab 100644 --- a/calvados/analysis.py +++ b/calvados/analysis.py @@ -28,9 +28,7 @@ def center_traj(pdb,traj,start=None,stop=None,step=1): """ Center trajectory """ - u = mda.Universe(pdb,traj) - with mda.Writer(f'{traj[:-4]}_c.dcd', len(u.atoms)) as W: for ts in u.trajectory[start:stop:step]: u.atoms.translate(-u.atoms.center_of_geometry() + 0.5 * u.dimensions[:3]) @@ -965,5 +963,6 @@ def calc_contact_map(path,sysname,output_path,chainid_dict={},is_slab=False,inpu cmap += (.5-.5*np.tanh((d-1.)/.3)).reshape(mask_frames.sum(), N_res_1,-1,N_res_2).sum(axis=0).sum(axis=1) cmap /= traj.n_frames + # save energy and contact maps np.save(output_path+f'/{sysname:s}_{name_1:s}_{name_2:s}_cmap.npy',cmap) diff --git a/calvados/build.py b/calvados/build.py index 2f9311b..10e21e0 100644 --- a/calvados/build.py +++ b/calvados/build.py @@ -1,7 +1,7 @@ from scipy import constants import numpy as np import pandas as pd -from openmm import unit +from openmm import app, unit from MDAnalysis import Universe from MDAnalysis.analysis import distances @@ -223,17 +223,17 @@ def build_xyzgrid(N,box): """ 3D grid """ r = box / np.sum(box) - a = np.cbrt(N / np.product(r)) + a = np.cbrt(N / np.prod(r)) n = a * r nxyz = np.floor(n) - while np.product(nxyz) < N: + while np.prod(nxyz) < N: ndeviation = n / nxyz devmax = np.argmax(ndeviation) nxyz[devmax] += 1 - while np.product(nxyz) > N: + while np.prod(nxyz) > N: nmax = np.argmax(nxyz) nxyz[nmax] -= 1 - if np.product(nxyz) < N: + if np.prod(nxyz) < N: nxyz[nmax] += 1 break @@ -298,7 +298,11 @@ def geometry_from_pdb(pdb,use_com=False): """ positions in nm""" with catch_warnings(): simplefilter("ignore") - u = Universe(pdb) + if pdb.lower().endswith('.cif'): + pdbx = app.pdbxfile.PDBxFile(pdb) + u = Universe(pdbx) + else: + u = Universe(pdb) ag = u.atoms ag.translate(-ag.center_of_mass()) if use_com: diff --git a/calvados/components.py b/calvados/components.py index 94f4362..af6ad58 100644 --- a/calvados/components.py +++ b/calvados/components.py @@ -8,6 +8,8 @@ import numpy as np +import os + class Component: """ Generic component of the system. """ @@ -33,7 +35,12 @@ def calc_comp_seq(self): """ Calculate sequence of component. """ if self.restraint: - self.seq, self.n_termini, self.c_termini = seq_from_pdb(f'{self.pdb_folder}/{self.name}.pdb') + input_pdb = f'{self.pdb_folder}/{self.name}.pdb' + input_pdbx = f'{self.pdb_folder}/{self.name}.cif' + if os.path.isfile(input_pdbx): + self.seq, self.n_termini, self.c_termini = seq_from_pdb(input_pdbx) + else: + self.seq, self.n_termini, self.c_termini = seq_from_pdb(input_pdb) else: records = read_fasta(self.ffasta) self.seq = str(records[self.name].seq) @@ -116,7 +123,12 @@ def calc_x_from_pdb(self): """ Calculate protein positions from pdb. """ input_pdb = f'{self.pdb_folder}/{self.name}.pdb' - self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdb + input_pdbx = f'{self.pdb_folder}/{self.name}.cif' + + if os.path.isfile(input_pdbx): + self.xinit, self.dimensions = build.geometry_from_pdb(input_pdbx,use_com=self.use_com) # read from pdb + else: + self.xinit, self.dimensions = build.geometry_from_pdb(input_pdb,use_com=self.use_com) # read from pdbx def calc_ssdomains(self): """ Get bounds for restraints (harmonic). """ diff --git a/calvados/sequence.py b/calvados/sequence.py index dc1c389..80aeffa 100644 --- a/calvados/sequence.py +++ b/calvados/sequence.py @@ -26,6 +26,8 @@ from calvados import analysis, interactions +from openmm import app + ### SEQUENCE INPUT / OUTPUT def read_fasta(ffasta): records = SeqIO.to_dict(SeqIO.parse(ffasta, "fasta")) @@ -36,7 +38,11 @@ def seq_from_pdb(pdb,selection='all',fmt='string'): with warnings.catch_warnings(): warnings.simplefilter("ignore") - u = Universe(pdb) + if pdb.lower().endswith('.cif'): + pdbx = app.pdbxfile.PDBxFile(pdb) + u = Universe(pdbx) + else: + u = Universe(pdb) # we do not assume residues in the PDB file are numbered from 1 n_termini = [0] diff --git a/calvados/sim.py b/calvados/sim.py index bebc9ef..110d773 100644 --- a/calvados/sim.py +++ b/calvados/sim.py @@ -215,9 +215,11 @@ def build_system(self): self.add_custom_restraints() self.pdb_cg = f'{self.path}/top.pdb' + self.cif_cg = f'{self.path}/top.cif' a = md.Trajectory(self.pos, self.top, 0, self.box, [90,90,90]) if self.restart != 'pdb': # only save new topology if no system pdb is given a.save_pdb(self.pdb_cg) + a.save_cif(self.cif_cg) self.add_forces_to_system() self.print_system_summary() @@ -507,10 +509,23 @@ def simulate(self): fcheck_out = f'{self.path}/restart.chk' append = False - if self.restart == 'pdb' and os.path.isfile(fcheck_in): - pdb = app.pdbfile.PDBFile(fcheck_in) + if self.restart == 'pdb': + if os.path.isfile(fcheck_in): + pdb = app.pdbfile.PDBFile(fcheck_in) + else: + pdb = app.pdbfile.PDBFile(self.pdb_cg) + elif self.restart == 'cif': + if os.path.isfile(fcheck_in): + pdb = app.pdbxfile.PDBxFile(fcheck_in) + else: + pdb = app.pdbxfile.PDBxFile(self.cif_cg) + elif self.restart == 'checkpoint': + if os.path.isfile(fcheck_in): + pdb = app.pdbxfile.PDBxFile(fcheck_in) + else: + pdb = app.pdbxfile.PDBxFile(self.cif_cg) else: - pdb = app.pdbfile.PDBFile(self.pdb_cg) + pdb = app.pdbxfile.PDBxFile(self.cif_cg) # use langevin integrator integrator = openmm.openmm.LangevinMiddleIntegrator(self.temp*unit.kelvin,self.friction_coeff/unit.picosecond,0.01*unit.picosecond) @@ -537,7 +552,7 @@ def simulate(self): print(f'Appending log file to {self.path}/{self.sysname:s}.log') simulation.loadCheckpoint(fcheck_in) else: - if self.restart == 'pdb': + if self.restart in ['pdb','cif']: print(f'Reading in system configuration {self.frestart}') elif self.restart == 'checkpoint': print(f'No checkpoint file {self.frestart} found: Starting from new system configuration') @@ -561,9 +576,11 @@ def simulate(self): simulation.reporters.append(app.dcdreporter.DCDReporter(f'{self.path}/equilibration_{self.sysname:s}.dcd',self.wfreq,append=append)) simulation.step(self.steps_eq) state_final = simulation.context.getState(getPositions=True) - rep = app.pdbreporter.PDBReporter(f'{self.path}/equilibration_final.pdb',0) - rep.report(simulation,state_final) - pdb = app.pdbfile.PDBFile(f'{self.path}/equilibration_final.pdb') + with open(f'{self.path}/equilibration_final.pdb', 'w') as f: + app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f) + with open(f'{self.path}/equilibration_final.cif', 'w') as f: + app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f) + pdb = app.pdbxfile.PDBxFile(f'{self.path}/equilibration_final.cif') for index, force in enumerate(self.system.getForces()): if isinstance(force, openmm.CustomExternalForce): @@ -584,9 +601,11 @@ def simulate(self): simulation.reporters.append(app.dcdreporter.DCDReporter(f'{self.path}/equilibration_{self.sysname:s}.dcd',self.wfreq,append=append)) simulation.step(self.steps_eq) state_final = simulation.context.getState(getPositions=True,enforcePeriodicBox=True) - rep = app.pdbreporter.PDBReporter(f'{self.path}/equilibration_final.pdb',0) - rep.report(simulation,state_final) - pdb = app.pdbfile.PDBFile(f'{self.path}/equilibration_final.pdb') + with open(f'{self.path}/equilibration_final.pdb', 'w') as f: + app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f) + with open(f'{self.path}/equilibration_final.cif', 'w') as f: + app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f) + pdb = app.pdbxfile.PDBxFile(f'{self.path}/equilibration_final.cif') topology = pdb.getTopology() a, b, c = state_final.getPeriodicBoxVectors() topology.setPeriodicBoxVectors(state_final.getPeriodicBoxVectors()) @@ -632,10 +651,10 @@ def simulate(self): dt_string = now.strftime("%Y%d%m_%Hh%Mm%Ss") state_final = simulation.context.getState(getPositions=True,enforcePeriodicBox=True) - rep = app.pdbreporter.PDBReporter(f'{self.path}/{self.sysname}_{dt_string}.pdb',0) - rep.report(simulation,state_final) - rep = app.pdbreporter.PDBReporter(f'{self.path}/checkpoint.pdb',0) - rep.report(simulation,state_final) + with open(f'{self.path}/checkpoint.pdb', 'w') as f: + app.PDBFile.writeFile(simulation.topology, state_final.getPositions(), f) + with open(f'{self.path}/checkpoint.cif', 'w') as f: + app.PDBxFile.writeFile(simulation.topology, state_final.getPositions(), f) def run(path='.',fconfig='config.yaml',fcomponents='components.yaml'): with open(f'{path}/{fconfig}','r') as stream: diff --git a/examples/slab_IDR/README.md b/examples/slab_IDR/README.md index a867a5b..e5b9fd2 100644 --- a/examples/slab_IDR/README.md +++ b/examples/slab_IDR/README.md @@ -1,7 +1,7 @@ The lines below run direct coexistence simulations of 100 copies of a single IDR: ```bash -python prepare.py --name --replica +python prepare.py --name --replica --gpu_id python _/run.py --path _ ``` diff --git a/examples/slab_IDR/prepare.py b/examples/slab_IDR/prepare.py index 68d395a..5ee71cf 100644 --- a/examples/slab_IDR/prepare.py +++ b/examples/slab_IDR/prepare.py @@ -34,7 +34,7 @@ wfreq = N_save, # dcd writing frequency, 1 = 10fs steps = N_frames*N_save, # number of simulation steps runtime = 0, # overwrites 'steps' keyword if > 0 - platform = 'CUDA', # 'CUDA' + platform = 'CPU', # 'CUDA' restart = 'checkpoint', frestart = 'restart.chk', verbose = True, diff --git a/examples/slab_IDR_MDP/prepare.py b/examples/slab_IDR_MDP/prepare.py index 658dd6e..bd9fac1 100644 --- a/examples/slab_IDR_MDP/prepare.py +++ b/examples/slab_IDR_MDP/prepare.py @@ -68,6 +68,10 @@ # calculate contact map between name_1 and name_1 chainid_dict = dict({args.name_1:s} = (0,99)) calc_contact_map(path="{path:s}",sysname="{sysname:s}",output_path="data",chainid_dict=chainid_dict,is_slab=True) + +# calculate contact map between name_2 and name_2 (note that this is not always meaningful) +# chainid_dict = dict({args.name_2:s} = (100,199)) +# calc_contact_map(path="{path:s}",sysname="{sysname:s}",output_path="data",chainid_dict=chainid_dict,is_slab=True) """ config.write(path,name='config.yaml',analyses=analyses) diff --git a/setup.py b/setup.py index 1adde58..1cc033c 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ setup( name='calvados', - version='0.8.0', + version='0.8.1', description='Coarse-grained implicit-solvent simulations of biomolecules', url='https://github.com/KULL-Centre/CALVADOS', authors=[ @@ -23,22 +23,23 @@ license='GNU GPL3', packages=find_packages(), install_requires=[ - 'OpenMM==8.2.0', - 'numpy==1.24', - 'pandas==2.1.1', - 'MDAnalysis==2.6.1', - 'biopython==1.81', - 'Jinja2==3.1.2', - 'tqdm==4.66', - 'matplotlib==3.8', - 'PyYAML==6.0', - 'statsmodels==0.14', - 'localcider==0.1.21', - 'pytest==8.3.3', - 'numba==0.60', - 'scipy==1.13', - 'mdtraj==1.10', + "OpenMM>=8.4,<8.5", + "MDAnalysis>=2.10,<2.11", + "mdtraj>=1.11,<1.12", + 'numpy', + 'pandas', + 'biopython', + 'Jinja2', + 'tqdm', + 'matplotlib', + 'PyYAML', + 'statsmodels', + 'localcider', + 'pytest', + 'numba', + 'scipy' ], + python_requires=">=3.13,<3.14", # include_package_data=True, package_data={'' : ['data/*.csv', 'data/*.yaml', 'data/templates/*']}, @@ -46,6 +47,6 @@ classifiers=[ 'Intended Audience :: Science/Research', 'Operating System :: POSIX :: Linux', - 'Programming Language :: Python :: 3>=3.7,<3.11', + 'Programming Language :: Python :: 3.13', ], )