Lincs warning during WT Metadynamics

User discussions
1

GROMACS version: 2021.6-Colvars-2022-12-16-dev
GROMACS modification: No
Hi all,
I am struggling with a simulation of a membrane protein with a small molecule: I generated the topology using CHARMM-GUI with AMBER ff’s (FF19SB, LIPID21 and GAFF2); the system has been equilibrated through several steps with decreasing restraints on different components (membrane and solvent, heavy atoms of protein and ligand) and then simulated without restraints in NPT. I then set up a well tempered metadynamics in which I have biased a distance and a dihedral angle with gaussian height varying from 6 to 1 kJ/mol and bias temperature from 6000 to 4000 in different attempts, but the simulation, after some nanoseconds, crashes with several lincs warnings; as far as I can tell the first warning always involves one of the atoms that define the biased dihedral. Below the colvars file I used in one of the attempts:

colvar {
name d
distance {
# heavy atoms of ligand
group1 { atomNumbers 19467 19470 19473 19474 19476 19477 19478 19482 19483 19484 19488 19490 19491 19494 19496 19499 19502 19503 19504 19507 19509 19511 19512 19513 19514 19515 19516 19517 19519 19521 19522 19524 19526 19527 19528 19532 19533 19534 19538 19539 19540
}
# C-alpha of THR46
group2 {atomNumbers 729
}
}
upperBoundary 3.2
width 0.05
}

colvar {
name tors_lig_rec
# Torsional defining orietation between central benzene ring of ligand and axis between ser51-thr46 C-alpha
dihedral {
group1 { atomNumbers 19499 }
group2 { atomNumbers 19512 }
group3 { atomNumbers 729 }
group4 { atomNumbers 796 }
}
}

metadynamics {
name meta_d_tors
colvars d tors_lig_rec
hillWeight 4
gaussianSigmas 0.05 5
outputEnergy on
wellTempered on
biasTemperature 6000
keepFreeEnergyFiles on
writeHillsTrajectory on
}
harmonicWalls {
name wall_d
colvars d
upperWalls 3.0
upperWallConstant 0.05
}

colvarsRestartFrequency 1000
colvarsTrajFrequency 1000

As an addendum regarding topology, I noticed that in the system generated by CHARMM-GUI there was a small net charge (something like -0.0001) coming from ligand topology, so I also tried to replace the topology of the ligand with one generated through acpype which gave a neutral system: equilibrations and plain md were successful but the issue with WT-MetaD is still there.

Does anyone have any hints about a possible solution?

Thanks
Stefano

2

Hi all,
I would like to add some details related to system instability described in the previous post: I tried to generate a “finer” parameterization of the ligand, through a reported optimization procedure based on energies and forces of about 10000 conformations obtained from dftb+ md, but I still got the same issue.
Without any bias I could simulate the system for up to 600 ns without any errors; moreover, when biasing only distance or a distance and an angle cv’s (involving two atoms of the ligand and one of the protein) the WT-Metad simulation runs smoothly.
Could the issue possibly be related to the definition of dihedral cv?

Thanks for any hint and advice
Stefano

3

It seems unlikely to be the small net charge.

I’m not familiar with colvars, but one troubleshooting approach could be to progressively, linearly ramp up the dihedral angle to know where the CV starts to be a problem, and/or whether it might instead be linked to the metaD parameters. (maybe a few replica for reproducibility)