search

bioexcel / biobb_wf_amber

This tutorials aim to illustrate the process of setting up a simulation system containing a protein, step by step, using the BioExcel Building Blocks library (biobb) wrapping the Ambertools MD package.
https://mmb.irbbarcelona.org/biobb/
Apache License 2.0
1 stars 3 forks source link

issues : ### Step 4: Generate ligand topology (parameters). #2

Closed vas2201 closed 1 year ago

vas2201 commented 1 year ago

Hi there,

While generating the ligand topology (parameters), I am getting the following error for the RNA-drug complex (6XB7). Could you please advise on this matter ?

Thank for your help.

Regards Vas

################code in juputer notebook ################

Create Ligand system topology, STEP 4

Acpype_params_gmx: Generation of topologies for AMBER with ACPype

Import module

from biobb_chemistry.acpype.acpype_params_ac import acpype_params_ac

Create prop dict and inputs/outputs

output_acpype_inpcrd = ligandCode+'params.inpcrd' output_acpype_frcmod = ligandCode+'params.frcmod' output_acpype_lib = ligandCode+'params.lib' output_acpype_prmtop = ligandCode+'params.prmtop' output_acpype = ligandCode+'params'

prop = { 'basename' : output_acpype, 'charge' : mol_charge }

Create and launch bb

acpype_params_ac(input_path=output_babel_min, output_path_inpcrd=output_acpype_inpcrd, output_path_frcmod=output_acpype_frcmod, output_path_lib=output_acpype_lib, output_path_prmtop=output_acpype_prmtop, properties=prop) ################################################################## Output & error :

Create Ligand system topology, STEP 4

Acpype_params_gmx: Generation of topologies for AMBER with ACPype

Import module

from biobb_chemistry.acpype.acpype_params_ac import acpype_params_ac

​

Create prop dict and inputs/outputs

output_acpype_inpcrd = ligandCode+'params.inpcrd'

output_acpype_frcmod = ligandCode+'params.frcmod'

output_acpype_lib = ligandCode+'params.lib'

output_acpype_prmtop = ligandCode+'params.prmtop'

output_acpype = ligandCode+'params'

​

prop = {

'basename' : output_acpype,

'charge' : mol_charge

}

​

Create and launch bb

acpype_params_ac(input_path=output_babel_min, 

            output_path_inpcrd=output_acpype_inpcrd,

            output_path_frcmod=output_acpype_frcmod,

            output_path_lib=output_acpype_lib,

            output_path_prmtop=output_acpype_prmtop,

            properties=prop)

2023-03-14 01:42:28,795 [MainThread ] [INFO ] acpype -i /home/nmrbox/spenumutchu/workflow2023/biobb_wf_amber_md_setup/biobb_wf_amber_md_setup/notebooks/mdsetup_lig/UYS.H.min.mol2 -b UYSparams.HndseR -n 0

2023-03-14 01:42:30,156 [MainThread ] [INFO ] Exit code 19

2023-03-14 01:42:30,157 [MainThread ] [INFO ] ========================================================================== | ACPYPE: AnteChamber PYthon Parser interfacE v. 2022.6.6 (c) 2023 AWSdS |

==> ... charge set to 0 ==> Executing Antechamber... ERROR: ++++++++++start_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ERROR: Welcome to antechamber 22.0: molecular input file processor.

Info: The atom type is set to gaff2; the options available to the -at flag are gaff, gaff2, amber, bcc, and sybyl.

Info: Finished reading file (UYS.H.min.mol2); atoms read (42), bonds read (43). Info: Determining atomic numbers from atomic symbols which are case sensitive. Running: /home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/bin/bondtype -j full -i ANTECHAMBER_BOND_TYPE.AC0 -o ANTECHAMBER_BOND_TYPE.AC -f ac Info: Bond types are assigned for valence state (1) with penalty (1).

Running: /home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/bin/atomtype -i ANTECHAMBER_AC.AC0 -o ANTECHAMBER_AC.AC -p gaff2 Info: Total number of electrons: 171; net charge: 0 Info: The number of electrons is odd (171). Please check the total charge (-nc flag) and spin multiplicity (-m flag).

Running: /home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/bin/sqm -O -i sqm.in -o sqm.out /home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/bin/wrapped_progs/antechamber: Fatal Error! Cannot properly run "/home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/bin/sqm -O -i sqm.in -o sqm.out".

ERROR: ++++++++++end_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ERROR: Antechamber failed ERROR: ++++++++++start_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ERROR: Cannot open file (UYSparams.HndseR_bcc_gaff2.mol2) with mode (r). No such file or directory ERROR: ++++++++++end_quote+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ERROR: Parmchk failed ERROR: Tleap failed ==> Removing temporary files... ERROR: ACPYPE FAILED: [Errno 2] No such file or directory: 'UYSparams.HndseR_AC.prmtop' Traceback (most recent call last): File "/home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/lib/python3.10/site-packages/acpype/cli.py", line 140, in init_main molecule.createMolTopol() File "/home/nmrbox/spenumutchu/anaconda3/envs/biobb_AMBER_MDsetup_tutorials/lib/python3.10/site-packages/acpype/topol.py", line 1116, in createMolTopol self.topFileData = open(self.acTopFileName).readlines() FileNotFoundError: [Errno 2] No such file or directory: 'UYSparams.HndseR_AC.prmtop' Total time of execution: 1s

2023-03-14 01:42:30,184 [MainThread ] [INFO ] Removed temporary folder: UYSparams.HndseR.acpype

gbayarri commented 1 year ago

Hi, it seems that you have passed a wrong charge parameter:

Please check the total charge (-nc flag) and spin multiplicity (-m flag).

This parameter is in the input properties for acpype and it must be assigned at the beginning of the workflow (mol_charge, below the ligand code).

There are tools for calculating this charge (ie, marvin suite) but not in the BioBB's, so you should find it by yourself.

Let me know if that worked for you.

vas2201 commented 1 year ago

Thank you ! it worked out.

How can I speed up the calculation?. Looks like it is only running on one CPU (total :128 cpus).

Thank you.

Import module from biobb_amber.sander.sander_mdrun import sander_mdrun

Create prop dict and inputs/outputs

output_nvt_traj_path = 'sander.nvt.netcdf' output_nvt_rst_path = 'sander.nvt.rst' output_nvt_log_path = 'sander.nvt.log'

prop = { "simulation_type" : 'nvt', "mdin" : { 'nstlim' : 500, # Reducing the number of steps for the sake of time (1ps) 'ntr' : 1, # Overwritting restrain parameter 'restraintmask' : '\"!:WAT,Cl-,Na+ & !@H=\"', # Restraining solute heavy atoms 'restraint_wt' : 5.0 # With a force constant of 5 Kcal/mol*A2 } }

Create and launch bb

sander_mdrun(input_top_path=output_ions_top_path, input_crd_path=output_heat_rst_path, input_ref_path=output_heat_rst_path, output_traj_path=output_nvt_traj_path, output_rst_path=output_nvt_rst_path, output_log_path=output_nvt_log_path, properties=prop)

log file :

Import module

from biobb_amber.sander.sander_mdrun import sander_mdrun

​

Create prop dict and inputs/outputs

output_nvt_traj_path = 'sander.nvt.netcdf'

output_nvt_rst_path = 'sander.nvt.rst'

output_nvt_log_path = 'sander.nvt.log'

​

prop = {

"simulation_type" : 'nvt',

"mdin" : { 

    'nstlim' : 500, # Reducing the number of steps for the sake of time (1ps)

    'ntr' : 1,     # Overwritting restrain parameter

    'restraintmask' : '\"!:WAT,Cl-,Na+ & !@H=\"',      # Restraining solute heavy atoms

    'restraint_wt' : 5.0                              # With a force constant of 5 Kcal/mol*A2

}

}

​

Create and launch bb

sander_mdrun(input_top_path=output_ions_top_path,

        input_crd_path=output_heat_rst_path,

        input_ref_path=output_heat_rst_path,

        output_traj_path=output_nvt_traj_path,

        output_rst_path=output_nvt_rst_path,

        output_log_path=output_nvt_log_path,

        properties=prop)

2023-03-15 02:35:00,708 [MainThread ] [INFO ] Executing biobb_amber.sander.sander_mdrun Version: 3.8.0 2023-03-15 02:35:00,713 [MainThread ] [INFO ] Creating 188dedfb-ccac-4c48-b0cf-2b5b6b5430ad temporary folder 2023-03-15 02:35:00,719 [MainThread ] [INFO ] sander -O -i 188dedfb-ccac-4c48-b0cf-2b5b6b5430ad/sander.mdin -p structure.ions.parmtop -c sander.heat.rst -r sander.nvt.rst -o sander.nvt.log -x sander.nvt.netcdf -ref sander.heat.rst

2023-03-15 02:36:47,593 [MainThread ] [INFO ] Exit code 0

2023-03-15 02:36:47,605 [MainThread ] [INFO ] Removed: ['188dedfb-ccac-4c48-b0cf-2b5b6b5430ad', 'mdinfo']

gbayarri commented 1 year ago

Hi, I'm glad that it worked for you.

Regarding to the CPU's, the Ambertools installed via conda for this BioBB doesn't include sander.MPI (required for running in more than one CPU). But if you have Ambertools installed in the computer where you are running this workflow, you can check if it's installed. In case you have it, you can run the workflow using your own Ambertools installation instead of the one included in the conda environment. For doing that, just add the binary_path property to the properties object of sander_mdrun:

prop = {
  "simulation_type" : 'nvt',
  "mdin" : { 
    'nstlim' : 500, # Reducing the number of steps for the sake of time (1ps)
    'ntr' : 1,     # Overwritting restrain parameter
    'restraintmask' : '\"!:WAT,Cl-,Na+ & !@H=\"',      # Restraining solute heavy atoms
    'restraint_wt' : 5.0                              # With a force constant of 5 Kcal/mol*A2
  },
  'binary_path': /path/to/sander/in/your/computer
}

Additionally, if you have an Amber license in the computer where you are executing the WF, you can use the pmemd_mdrun module instead of the sander_mdrun in the same way (adding the binary_path property):

https://biobb-amber.readthedocs.io/en/latest/pmemd.html#module-pmemd.pmemd_mdrun

vas2201 commented 1 year ago

Thanks for your message. I update code as you suggested, However it is not using sander.MPI module. Still using sander module . I did not see any errors. I also confirmed from top command that 'sander ' is running.

if possible , please advise on this matter.

Import module

from biobb_amber.sander.sander_mdrun import sander_mdrun

Create prop dict and inputs/outputs

output_heat_traj_path = 'sander.heat.netcdf' output_heat_rst_path = 'sander.heat.rst' output_heat_log_path = 'sander.heat.log'

prop = { "simulation_type" : "heat", "mdin" : { 'nstlim' : 2500, # Reducing the number of steps for the sake of time (5ps) 'ntr' : 1, # Overwritting restrain parameter 'restraintmask' : '\"!:WAT,Cl-,Na+\"', # Restraining solute 'restraint_wt' : '10.0', # With a force constant of 10 Kcal/mol*A2 }, 'binary_path' : '/usr/software/amber/bin/sander.MPI' }

Create and launch bb

sander_mdrun(input_top_path=output_ions_top_path, input_crd_path=output_min_rst_path, input_ref_path=output_min_rst_path, output_traj_path=output_heat_traj_path, output_rst_path=output_heat_rst_path, output_log_path=output_heat_log_path, properties=prop)

gbayarri commented 1 year ago

Hi, can you please provide us the complete log for this step?

If you see the log you put in your yesterday's post, you will see this line:

2023-03-15 02:35:00,719 [MainThread ] [INFO ] sander -O -i 188dedfb-ccac-4c48-b0cf-2b5b6b5430ad/sander.mdin -p structure.ions.parmtop -c sander.heat.rst -r sander.nvt.rst -o sander.nvt.log -x sander.nvt.netcdf -ref sander.heat.rst

Now, with the binary_path, you should see the complete path from your computer instead of sander:

2023-03-15 02:35:00,719 [MainThread ] [INFO ] /usr/software/amber/bin/sander.MPI -O -i 188dedfb-ccac-4c48-b0cf-2b5b6b5430ad/sander.mdin -p structure.ions.parmtop -c sander.heat.rst -r sander.nvt.rst -o sander.nvt.log -x sander.nvt.netcdf -ref sander.heat.rst
vas2201 commented 1 year ago

Hi,

Thanks for your reply.

output from below cell of NVT

2023-03-16 04:55:08,532 [MainThread ] [INFO ] Executing biobb_amber.sander.sander_mdrun Version: 3.8.0 2023-03-16 04:55:08,538 [MainThread ] [INFO ] Creating 4ec74b60-6d82-4420-b9c1-214fcfd1e463 temporary folder 2023-03-16 04:55:08,545 [MainThread ] [INFO ] sander -O -i 4ec74b60-6d82-4420-b9c1-214fcfd1e463/sander.mdin -p structure.ions.parmtop -c sander.min.rst -r sander.heat.rst -o sander.heat.log -x sander.heat.netcdf -ref sander.min.rst

2023-03-16 04:56:49,372 [MainThread ] [INFO ] Exit code 0

2023-03-16 04:56:49,391 [MainThread ] [INFO ] Removed: ['4ec74b60-6d82-4420-b9c1-214fcfd1e463', 'mdinfo']

0

sander.nvt.log

      -------------------------------------------------------
      Amber 22 SANDER                              2022
      -------------------------------------------------------

| Run on 03/16/2023 at 04:56:49

| Executable path: sander | Working directory: /home/nmrbox/spenumutchu/workflow2023/biobb_wf_amber_md_setup/biobb_wf_amber_md_setup/notebooks/mdsetup/hnrnp_A1 | Hostname: indium.nmrbox.org

[-O]verwriting output

File Assignments: | MDIN: b9f9752f-410a-4ea4-a600-5b231b4073ff/sander.mdin
| MDOUT: sander.nvt.log
|INPCRD: sander.heat.rst
| PARM: structure.ions.parmtop
|RESTRT: sander.nvt.rst
| REFC: sander.heat.rst
| MDVEL: mdvel
| MDFRC: mdfrc
| MDEN: mden
| MDCRD: sander.nvt.netcdf
|MDINFO: mdinfo
| MTMD: mtmd
|INPDIP: inpdip
|RSTDIP: rstdip
|INPTRA: inptraj

Here is the input file:

This mdin file has been created by the biobb_amber module from the BioBB library Type of mdin: nvt
&cntrl
imin = 0 ! BioBB simulation_type nvt|npt|free|heat
cut = 10.0 ! BioBB simulation_type nvt|npt|free|heat
ntr = 1 ! BioBB property
ntc = 2 ! BioBB simulation_type nvt|npt|free|heat
ntf = 2 ! BioBB simulation_type nvt|npt|free|heat
ntt = 3 ! BioBB simulation_type nvt|npt|free|heat
ig = -1 ! BioBB simulation_type nvt|npt|free|heat
ioutfm = 1 ! BioBB simulation_type nvt|npt|free|heat
iwrap = 1 ! BioBB simulation_type nvt|npt|free|heat
nstlim = 500 ! BioBB property
dt = 0.002 ! BioBB simulation_type nvt|npt|free|heat
irest = 1 ! BioBB simulation_type nvt
gamma_ln = 5.0 ! BioBB simulation_type nvt
ntb = 1 ! BioBB simulation_type nvt
ntx = 5 ! BioBB simulation_type nvt
restraintmask = "!:WAT,Cl-,Na+ & !@H=" ! BioBB property
restraint_wt = 5.0 ! BioBB property
&end

  1. RESOURCE USE:

| Flags:
getting box info from netcdf restart file | NetCDF restart box info found |Largest sphere to fit in unit cell has radius = 55.790 | New format PARM file being parsed. | Version = 1.000 Date = 03/16/23 Time = 04:45:42 NATOM = 179870 NTYPES = 17 NBONH = 176810 MBONA = 2796 NTHETH = 5654 MTHETA = 3729 NPHIH = 12197 MPHIA = 11271 NHPARM = 0 NPARM = 0 NNB = 261143 NRES = 58801 NBONA = 2796 NTHETA = 3729 NPHIA = 11271 NUMBND = 67 NUMANG = 152 NPTRA = 192 NATYP = 37 NPHB = 0 IFBOX = 2 NMXRS = 24 IFCAP = 0 NEXTRA = 0 NCOPY = 0

| Memory Use Allocated | Real 13227071 | Hollerith 598413 | Integer 5221334 | Max Pairs 103605120 | nblistReal 2158440 | nblist Int 6011970 | Total 571124 kbytes

| Note: 1-4 EEL scale factors are being read from the topology file.

| Note: 1-4 VDW scale factors are being read from the topology file. | Duplicated 0 dihedrals | Duplicated 0 dihedrals

 BOX TYPE: TRUNCATED OCTAHEDRON

Note: ig = -1. Setting random seed to 847180 based on wallclock time in microseconds.

  1. CONTROL DATA FOR THE RUN

default_name

General flags: imin = 0, nmropt = 0

Nature and format of input: ntx = 5, irest = 1, ntrx = 1

Nature and format of output: ntxo = 2, ntpr = 50, ntrx = 1, ntwr = 500 iwrap = 1, ntwx = 0, ntwv = 0, ntwe = 0 ioutfm = 1, ntwprt = 0, idecomp = 0, rbornstat= 0

Potential function: ntf = 2, ntb = 1, igb = 0, nsnb = 25 ipol = 0, gbsa = 0, iesp = 0 dielc = 1.00000, cut = 10.00000, intdiel = 1.00000

Frozen or restrained atoms: ibelly = 0, ntr = 1 restraint_wt = 5.00000

Molecular dynamics: nstlim = 500, nscm = 0, nrespa = 1 t = 0.00000, dt = 0.00200, vlimit = 20.00000

Langevin dynamics temperature regulation: ig = 847180 temp0 = 300.00000, tempi = 0.00000, gamma_ln= 5.00000

SHAKE: ntc = 2, jfastw = 0 tol = 0.00001

Ewald parameters: verbose = 0, ew_type = 0, nbflag = 1, use_pme = 1 vdwmeth = 1, eedmeth = 1, netfrc = 1 Box X = 136.656 Box Y = 136.656 Box Z = 136.656 Alpha = 109.471 Beta = 109.471 Gamma = 109.471 NFFT1 = 144 NFFT2 = 144 NFFT3 = 144 Cutoff= 10.000 Tol =0.100E-04 Ewald Coefficient = 0.27511 Interpolation order = 4 | INFO: Old style inpcrd file read

LOADING THE CONSTRAINED ATOMS AS GROUPS

  1. REFERENCE ATOM COORDINATES

    defa Mask !:WAT,Cl-,Na+ & !@H=; matches 2735 atoms

  1. ATOMIC COORDINATES AND VELOCITIES

default_name
begin time read from input coords = 0.100 ps

Number of triangulated 3-point waters found: 58104

 Sum of charges from parm topology file =  -0.00000024
 Forcing neutrality...

  1. RESULTS

| # of SOLUTE degrees of freedom (RNDFP): 362800. | # of SOLVENT degrees of freedom (RNDFS): 0. | NDFMIN = 362800. NUM_NOSHAKE = 0 CORRECTED RNDFP = 362800. | TOTAL # of degrees of freedom (RNDF) = 362800.

APPROXIMATING switch and d/dx switch using CUBIC SPLINE INTERPOLATION using 5000.0 points per unit in tabled values TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff | CHECK switch(x): max rel err = 0.2738E-14 at 2.422500 | CHECK d/dx switch(x): max rel err = 0.8314E-11 at 2.736960

| Local SIZE OF NONBOND LIST = 59686796 | TOTAL SIZE OF NONBOND LIST = 59686796

NSTEP = 50 TIME(PS) = 0.200 TEMP(K) = 157.74 PRESS = 0.0 Etot = -595059.2001 EKtot = 56860.9861 EPtot = -651920.1861 BOND = 641.9514 ANGLE = 1712.8105 DIHED = 4319.1981 1-4 NB = 1074.8502 1-4 EEL = 14053.5950 VDWAALS = 104676.6294 EELEC = -778743.3018 EHBOND = 0.0000 RESTRAINT = 344.0810 EAMBER (non-restraint) = -652264.2671 Ewald error estimate: 0.8579E-04

NSTEP = 100 TIME(PS) = 0.300 TEMP(K) = 211.11 PRESS = 0.0 Etot = -553784.5931 EKtot = 76099.6110 EPtot = -629884.2041 BOND = 828.2652 ANGLE = 2093.4503 DIHED = 4345.9438 1-4 NB = 1106.6725 1-4 EEL = 14029.9697 VDWAALS = 101575.7710 EELEC = -754383.4902 EHBOND = 0.0000 RESTRAINT = 519.2137 EAMBER (non-restraint) = -630403.4178 Ewald error estimate: 0.6929E-04

NSTEP = 150 TIME(PS) = 0.400 TEMP(K) = 243.15 PRESS = 0.0 Etot = -527830.5039 EKtot = 87649.4355 EPtot = -615479.9394 BOND = 921.0015 ANGLE = 2410.1813 DIHED = 4477.4123 1-4 NB = 1120.5144 1-4 EEL = 13988.9915 VDWAALS = 96508.2033 EELEC = -735391.6873 EHBOND = 0.0000 RESTRAINT = 485.4436 EAMBER (non-restraint) = -615965.3830 Ewald error estimate: 0.1773E-04

log14.out

2023-03-16 04:56:49,848 [MainThread ] [INFO ] Executing biobb_amber.sander.sander_mdrun Version: 3.8.0 2023-03-16 04:56:49,854 [MainThread ] [INFO ] Creating b9f9752f-410a-4ea4-a600-5b231b4073ff temporary folder 2023-03-16 04:56:49,860 [MainThread ] [INFO ] sander -O -i b9f9752f-410a-4ea4-a600-5b231b4073ff/sander.mdin -p structure.ions.parmtop -c sander.heat.rst -r sander.nvt.rst -o sander.nvt.log -x sander.nvt.netcdf -ref sander.heat.rst

gbayarri commented 1 year ago

Hi, I see that you are running the version 3.8.0, the binary_path property was introduced in the 3.9.0 version. In your case, the property name should be sander_path instead of binary_path.

So, please try with this new property:

prop = {
  "simulation_type" : 'nvt',
  "mdin" : { 
    'nstlim' : 500, # Reducing the number of steps for the sake of time (1ps)
    'ntr' : 1,     # Overwritting restrain parameter
    'restraintmask' : '\"!:WAT,Cl-,Na+ & !@H=\"',      # Restraining solute heavy atoms
    'restraint_wt' : 5.0                              # With a force constant of 5 Kcal/mol*A2
  },
  'sander_path': /path/to/sander/in/your/computer
}

If this doesn't work either, add this propery and try again:

mpi_bin: mpirun

vas2201 commented 1 year ago

Thank you ! That worked for me. I do have another request for you. We usually use emap restraints to do the MD simulations. How can I implement the below code to biobb_wf_amber_md_setup.ipynab file.
Really appreciate all your help on this matter.

&cntrl nstlim=2500000, # number of MD steps iemap=1 # specify electrostatic map calculations cut=999.9, # non-bonded cutoff igb=1, # implicit solvent model rgbmax=10, # maximum value of solvent radius saltcon=0.155, # salt concentration ntx=1, # initial velocities generated using random number generator irest=0, # start a new simulation nscm=1000, # frequency for the SHAKE algorithm to converge tempi=300.0, # starting temperature temp0=300.0, # reference temperature ntt=3, # Nose-Hoover thermostat gamma_ln=2.0, # Nose-Hoover relaxation time ig=-1, # exclude atoms from IG list ntpr=5000, # frequency to write output to mdout file ntwx=5000, # frequency to write trajectory file ntwr=500000, # frequency to write restart file nmropt=1, # optimize hydrogen bond network dt=0.002, # time step size ntc=2, ntf=2, # SHAKE algorithm applied to both bonds and angles, no constraints applied to dihedrals ioutfm=1, # format of the mdout file /

&emap # electrostatic map options mapfile=damfilt.pdb, atmask=':*', fcons=0.005, # electrostatic map input options resolution=5, # grid resolution for the electrostatic map mapfit='', # no fitting for the electrostatic map molfit='', # no fitting for the molecule /

0 to 300K (100 ps) 300 to 300 (500 ps) 300 to 0 (100 ps) 0 to 0 (50 ps)

&wt type='TEMP0', istep1=0, istep2=2500000, value1=300.0, value2=300.0,/

Simple constant temp incorporation of restraints

&wt type='REST', istep1=0,istep2=2500000,value1=1.0,value2=1.0,/

&wt type ='END' # end of weight options /

LISTOUT=POUT # write output to POUT file DISANG=rnaRST # restraint file for distance restraints END # end of input file

amberprod_5ns.in # append the input to the amberprod_5ns.in file

run the simulation using pmemd.cuda program

$AMBERHOME/bin/pmemd.cuda -O -i amberprod_5ns.in -o ${SLURM_ARRAY_TASK_ID}sim.out -p ${SLURM_ARRAY_TASK_ID}.prmtop -c ${SLURM_ARRAY_TASK_ID}min.rst -r ${SLURM_ARRAY_TASK_ID}sim.rst -ref ${SLURM_ARRAY_TASK_ID}min.rst -x ${SLURM_ARRAY_TASK_ID}sim.trj

write output to specified files

cp -ru $SLURM_SUBMIT_DIR rm -rf "$PFSDIR"/

gbayarri commented 1 year ago

Hi, in the current version of biobb_amber, this is not possible.

There is an input_mdin_path where you can put the path of a .mdin file, but then we parse it and it only would work until the first slash. That is:

&cntrl
nstlim=2500000, # number of MD steps
iemap=1 # specify electrostatic map calculations
cut=999.9, # non-bonded cutoff
igb=1, # implicit solvent model
rgbmax=10, # maximum value of solvent radius
saltcon=0.155, # salt concentration
ntx=1, # initial velocities generated using random number generator
irest=0, # start a new simulation
nscm=1000, # frequency for the SHAKE algorithm to converge
tempi=300.0, # starting temperature
temp0=300.0, # reference temperature
ntt=3, # Nose-Hoover thermostat
gamma_ln=2.0, # Nose-Hoover relaxation time
ig=-1, # exclude atoms from IG list
ntpr=5000, # frequency to write output to mdout file
ntwx=5000, # frequency to write trajectory file
ntwr=500000, # frequency to write restart file
nmropt=1, # optimize hydrogen bond network
dt=0.002, # time step size
ntc=2, ntf=2, # SHAKE algorithm applied to both bonds and angles, no constraints applied to dihedrals
ioutfm=1, # format of the mdout file
/

The rest of the parameters won't be read.

We opened an issue for adding a new property in the next release in order to avoid this parsing, and leave the .mdin file as it is:

https://github.com/bioexcel/biobb_amber/issues/6

gbayarri commented 1 year ago

Hi, it has been fixed in biobb_amber 4.0.0, recently released.

vas2201 commented 1 year ago

Thanks for your help. If possible, Could you please help me to provide an example to work with emap restraints ?. Regards Vas