marrink-lab / vermouth-martinize

Describe and apply transformation on molecular structures and topologies
Apache License 2.0
86 stars 38 forks source link

Potential problem with "-el" flag #479

Closed MikkelDA closed 1 year ago

MikkelDA commented 1 year ago

These problems were encountered in the latest development version of martinize2/vermouth that was available on the 18/10-2022.

I think i have encountered a bug with the "-el" flag or maybe my understanding of it is simply incorrect. My understanding is that the "-el" flag sets the lower limit for how long a Rubber band can be (while "-eu" sets the upper limit), meaning that the following must be true for "Rubber bands": el < Bond length < eu However i just noticed that some of the networks generated for my structure have bond lengths less than "-el", while none of them are above "-eu". Specifically 290 out of 2428 rubber bands had lengths less than "-el"

This is my martinize2 command: martinize2 -f 03_prot_editconf.pdb -x 04_output_martinize.pdb -o topol_res152_standard.top -ff martini3001 -from charmm -dssp /mnt/software/dssp/3.1.2/bin/mkdssp -p backbone -elastic -ef 700 -el 0.5 -eu 0.9 -ea 0 -ep 0 -scfix -cys auto -resid input -mutate HSE:HSD

The shortest distance martinize2 generated was: 0.38086 (Note that the longest distance created is 0.89991, indicating that "eu" is a hard cutoff, which i assumed "el" to also be).

Is this working as intended? Am i understanding it wrong?

pckroon commented 1 year ago

Thanks for reaching out. I think the issue here is in your interpretation (and thereby insufficient documentation from our side). From martinize2 -h: -el RB_LOWER_BOUND Elastic bond lower cutoff: F = Fc if rij < lo (default: 0) So, for all elastic bonds shorter than 0.5 in your example the force constant will be 700.

For all elastic bonds between -el and -eu (0.5 and 0.9) the force constant will be modulated with np.exp(-rate * ((distance - shift) ** power)) (https://github.com/marrink-lab/vermouth-martinize/blob/master/vermouth/processors/apply_rubber_band.py#L61). Rate, shift and power are decay_factor (-ea), lower_bound (-el), and decay_power respectively (-ep),

MikkelDA commented 1 year ago

Thanks for the reply