GROMACS version: 2023.2
GROMACS modification: No
Protein FF: Amber14sb
Ligand FF: GAFF2
Water FF: TIP3P
Simulation Type: Relative binding free energy calculations using Non Equilibrium TI method
I am trying to run an alchemical calculation to estimate the relative binding free energy difference between two ligands for a protein receptor. The protocol I am following is the one mentioned in (Gapsys et. al, Large scale relative protein ligand binding affinities using non-equilibrium alchemy). I am using the Amber14sb forcefield for the protein and GAFF2 for the ligands. The protocol first employs energy minimization on the system containing the protein with the ligand with dummy atoms for both states A and B. The last frame from minimization is then used to setup a short NVT run for equilibration in both states. The published protocol places Lincs constraints on all bonds for these dynamics simulations. I get multiple warnings here starting from the first frame itself. Some of the errors are listed below,
Errors
Command line:
gmx mdrun -ntomp 32 -ntmpi 1 -s tpr.tpr
Reading file tpr.tpr, VERSION 2023.2 (single precision)
Changing nstlist from 10 to 100, rlist from 1.1 to 1.262
Update groups can not be used for this system because there are three or more consecutively coupled constraints
1 GPU selected for this run.
Mapping of GPU IDs to the 2 GPU tasks in the 1 rank on this node:
PP:0,PME:0
PP tasks will do (non-perturbed) short-ranged interactions on the GPU
PP task will update and constrain coordinates on the CPU
PME tasks will do all aspects on the GPU
Using 1 MPI thread
Using 32 OpenMP threads
starting mdrun ‘protein and ligand in water’
3000000 steps, 6000.0 ps.
Step 0, time 0 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000758, max 0.040537 (between atoms 4696 and 4702)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
4696 4702 39.7 0.1797 0.1455 0.1399
Step 0, time 0 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000614, max 0.033828 (between atoms 4696 and 4702)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
4696 4702 40.8 0.1797 0.1446 0.1399
Step 0 Warning: pressure scaling more than 1%, mu: 0.989131 0.989131 0.989131
Step 1, time 0.002 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1246.896606, max 41854.613281 (between atoms 372 and 373)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
347 348 31.2 0.1078 0.1312 0.1090
347 360 35.6 0.1506 0.1995 0.1522
360 361 47.6 0.1216 0.1854 0.1229
360 362 83.1 0.1321 0.2714 0.1335
362 363 85.6 0.0999 0.1881 0.1010
362 364 94.9 0.1434 0.6319 0.1449
364 365 96.9 0.1078 0.2964 0.1090
364 366 94.7 0.1510 0.6541 0.1526
364 380 92.8 0.1506 0.5101 0.1522
366 367 94.8 0.1078 0.6398 0.1090
366 368 94.2 0.1079 0.8516 0.1090
366 369 78.3 0.1494 0.8606 0.1510
369 370 127.0 0.1386 924.4766 0.1400
369 378 90.5 0.1385 1.1291 0.1400
370 371 126.9 0.1068 923.5626 0.1080
372 373 132.7 0.1068 4520.4058 0.1080
372 374 113.6 0.1385 4521.9683 0.1400
374 375 101.1 0.1068 0.2392 0.1080
374 376 55.5 0.1385 0.1307 0.1400
376 377 58.8 0.1068 0.1775 0.1080
376 378 94.9 0.1385 0.3829 0.1400
378 379 94.5 0.1069 0.3353 0.1080
380 381 78.2 0.1216 0.1646 0.1229
380 382 80.0 0.1321 0.2305 0.1335
382 383 54.7 0.0999 0.1707 0.1010
382 384 51.0 0.1433 0.2323 0.1449
384 385 33.2 0.1078 0.1336 0.1090
384 386 30.7 0.1078 0.1304 0.1090
4692 4693 48.5 0.1522 0.2284 0.1538
4693 4694 87.0 0.1451 0.2804 0.1466
4693 4718 56.1 0.1087 0.1879 0.1098
4693 4719 56.5 0.1087 0.1899 0.1098
4694 4695 93.8 0.1454 0.6055 0.1463
4694 4703 85.0 0.1451 0.2794 0.1466
4695 4720 97.8 0.1091 0.2837 0.1097
4695 4721 97.9 0.1091 0.3791 0.1097
4696 4697 107.8 0.1387 0.7706 0.1399
4696 4702 119.1 0.1430 4159.7456 0.1399
4697 4698 106.0 0.1387 0.3090 0.1399
4697 4722 53.9 0.1077 0.0995 0.1086
4698 4699 98.1 0.1388 1.3131 0.1399
4698 4723 111.5 0.1076 0.2818 0.1086
4699 4700 89.9 0.1388 2.5366 0.1399
4699 4724 99.7 0.1078 0.8065 0.1086
4700 4701 100.6 0.1345 1.0309 0.1350
4700 4702 56.4 0.1408 4164.6294 0.1399
4702 4725 134.9 0.1080 4156.9385 0.1086
4703 4704 57.0 0.1522 0.2614 0.1538
4703 4726 69.8 0.1087 0.2384 0.1098
4703 4727 64.9 0.1086 0.2213 0.1098
4704 4728 33.5 0.1086 0.1328 0.1098
4704 4729 34.9 0.1086 0.1344 0.1098
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
4700 4734 124.4 0.1085 1.2337 0.1086
Back Off! I just backed up step1b.pdb to ./#step1b.pdb.6#
Back Off! I just backed up step1c.pdb to ./#step1c.pdb.6#
Wrote pdb files with previous and current coordinates
Step 1, time 0.002 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1246.572266, max 41847.437500 (between atoms 372 and 373)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
345 347 39.9 0.1433 0.1843 0.1449
347 348 45.9 0.1078 0.1487 0.1090
347 349 37.8 0.1509 0.1895 0.1526
347 360 44.8 0.1506 0.2057 0.1522
360 361 64.1 0.1216 0.2040 0.1229
360 362 93.6 0.1321 0.3628 0.1335
362 363 94.6 0.0999 0.2385 0.1010
362 364 98.2 0.1434 0.8354 0.1449
364 365 105.3 0.1078 0.3466 0.1090
364 366 102.9 0.1510 0.7213 0.1526
364 380 100.4 0.1506 0.6504 0.1522
366 367 99.5 0.1078 0.7380 0.1090
366 368 96.6 0.1079 0.9864 0.1090
366 369 96.9 0.1494 0.7140 0.1510
369 370 127.0 0.1386 924.7332 0.1400
369 378 99.4 0.1385 0.8317 0.1400
370 371 126.9 0.1068 923.8002 0.1080
372 373 132.8 0.1068 4519.6318 0.1080
372 374 113.6 0.1385 4520.8784 0.1400
374 375 109.5 0.1068 0.2936 0.1080
374 376 122.7 0.1385 0.2299 0.1400
376 377 64.9 0.1068 0.1426 0.1080
376 378 110.9 0.1385 0.3820 0.1400
378 379 106.0 0.1069 0.3559 0.1080
380 381 87.0 0.1216 0.1995 0.1229
380 382 91.9 0.1321 0.3400 0.1335
382 383 72.3 0.0999 0.1812 0.1010
382 384 73.1 0.1433 0.2819 0.1449
384 385 53.3 0.1078 0.1627 0.1090
384 386 52.9 0.1078 0.1650 0.1090
384 387 46.3 0.1506 0.2089 0.1522
4675 4676 37.4 0.1526 0.2010 0.1543
4676 4677 35.4 0.1452 0.1920 0.1468
4676 4692 50.6 0.1522 0.2202 0.1538
4676 4704 88.3 0.1522 0.2973 0.1538
4692 4693 99.6 0.1522 0.3684 0.1538
4692 4716 67.8 0.1086 0.1883 0.1098
4692 4717 71.0 0.1086 0.1942 0.1098
4693 4694 92.6 0.1451 0.9248 0.1466
4693 4718 105.1 0.1087 0.3343 0.1098
4693 4719 104.5 0.1087 0.3620 0.1098
4694 4695 91.4 0.1454 2.1200 0.1463
4694 4703 87.4 0.1451 0.9157 0.1466
4695 4696 50.2 0.1503 1.1733 0.1515
4695 4720 92.9 0.1091 1.1450 0.1097
4695 4721 93.0 0.1091 1.4526 0.1097
4696 4697 95.7 0.1387 2.5934 0.1399
4696 4702 119.2 0.1430 4159.3696 0.1399
4697 4698 97.9 0.1387 0.3180 0.1399
4697 4722 88.4 0.1077 0.1212 0.1086
4698 4699 97.5 0.1388 1.2543 0.1399
4698 4723 113.7 0.1076 0.2758 0.1086
4699 4700 90.6 0.1388 3.3103 0.1399
4699 4724 98.9 0.1078 0.7612 0.1086
4700 4701 97.7 0.1345 0.9151 0.1350
4700 4702 56.4 0.1408 4162.8477 0.1399
4702 4725 134.9 0.1080 4154.6221 0.1086
4703 4704 104.3 0.1522 0.4865 0.1538
4703 4726 99.9 0.1087 0.5108 0.1098
4703 4727 100.4 0.1086 0.4395 0.1098
4704 4728 91.8 0.1086 0.2760 0.1098
4704 4729 92.1 0.1086 0.2728 0.1098
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
4700 4734 118.2 0.1085 1.0937 0.1086
Wrote pdb files with previous and current coordinates
WARNING: Listed nonbonded interaction between particles 4693 and 4696
at distance 3.150 which is larger than the table limit 2.262 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
WARNING: Listed nonbonded interaction between particles 364 and 370
at distance 924.579 which is larger than the table limit 2.262 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
Program: gmx mdrun, version 2023.2
Source file: src/gromacs/gpu_utils/cudautils.cuh (line 190)
Fatal error:
Unexpected cudaStreamQuery failure. CUDA error #700 (cudaErrorIllegalAddress):
an illegal memory access was encountered.
For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage https://www.gromacs.org documentation
The NVT simulation then fails. If I only constrain the hydrogen bonds for these runs I don’t get these errors and the NVT equilibration and the NPT production run executes without errors. The non equilibrium transition runs also execute without errors but I can only constrain hydrogen bonds in all of these simulations. At any step if I try to constrain all bonds (which is what the protocol in the paper did) I get the same errors as posted above. Also, I get these errors only for the protein-ligand complex leg of the simulations and the ligand-water leg simulations run fine with constraints on all bonds. The protein structure itself is also stable and I have run serial MD simulations (following Prof. Lemkul’s protocols, constraints only on H-bonds) of the protein receptor in complex with the individual ligand poses (which seem to be stable) with the same forcefield combination without errors (overall RMSD of ~2 Angstroms over 200ns).
I am attaching the energy minimization mdp file and the NVT equilibration mdp file (for state A) here if it helps catch any mistakes that I am making in the input parameters. If it helps, I can also share the structure files for the energy minimization and NVT equilibration runs but for now I would prefer that I keep that under wraps.
I understand that these questions regarding Lincs warning have been asked before and I apologize for repeating these queries here but I wanted to confirm that the setup I am using is sensible. These calculations are quite expensive and involved and I would like to ensure that I am not wasting compute and time with incorrect simulation criteria.
I wanted to ask if it is feasible/recommended to run RBFE calculations with constraints only on h-bonds? For some other ligand pairs with the same receptor structure I can constrain all bonds and still get all steps of the protocol to run so is this an issue with the initial binding pose that I am using for the ligands and I should refine them further? It would be very helpful if I could get some clarity on this. Thank you!
nvt.mdp (2.6 KB)
em_l0.mdp (9.9 KB)