From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Thu Feb 14 2019 - 09:59:17 CST
Hi Hristina,
There are two potentially problematic areas:
1. Sampling time. Each sample is being run for 1ps based on your inputs. I don't know where your references are coming from, but especially for small systems, sampling times that are 100x longer would be where I start to trust the first or maybe the second significant digit of the result.
2. Probably less of a big deal here since its just in water, but when you are looking at ion exchanges in your protein, the end points are going to be very slow to converge, since the "disappeared" ion will be free to float out of the protein. Usually, in dual topology systems, the shared atoms prevent this from happening, and are why patches rather than full-blown new residues are the preferred approach. For your systems, consider pinning together the appearing/disappearing ions using extrabonds or colvars. There is a bit more reading in the mailing list archive that might be of interest: https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2017-2018/0465.html
-Josh
On 2019-02-13 18:03:26-07:00 owner-namd-l_at_ks.uiuc.edu wrote:
Hello everyone,
I have the following issue. A protein I'm interested in prefers to bind some anions over other anions. I want to evaluate the preference of the binding site for the different anions. The anions may have different charges and numbers of atoms, e.g. Cl- to HCO3- or HCO3- to CO3^2-. I tried to follow the NAMD tutorial for alchemical perturbations (the part for mutation of Tyr to Ala) with dual topology for these two alchemical transformations in a water box (no protein, only ion in about 3000 water molecules) to see if I'm getting sensible numbers. I have free energies of hydration for all 3 anions from previous calculations that were done with CHARMM, so I was comparing the DeltaDeltaG(hydration) values with the DeltaDeltaG(alchemical perturbation) values. Presumably, they should be the same or very similar (from thermodynamic cycle considerations). However, I'm getting significant difference (~10 kcal/mol), which means that I'm doing something wrong. I've used the same force fields parameters in all calculations (all parameters are available in the current CHARMM toppar files). The input files were modified from the tutorial files or Tyr to Ala alchemical perturbation. I haven't changed anything in the FEP block. These are the dual topologies I used:
RESI ZERO -1.00
GROUP
ATOM CA CG2O6 0.69 ! H3 O1
ATOM H3A HGP1 0.43 ! | /
ATOM O1A OG2D2 -0.76 ! O3--C (-)
ATOM O2A OG2D2 -0.76 ! Cl \
ATOM O3A OG311 -0.60 ! O2
GROUP
ATOM CLB CLA -1.00
BOND CA O1A CA O2A CA O3A O3A H3A
IMPR CA O1A O2A O3A
! seed = O1A CA O2A
IC O2A O1 *CA O3A 0.0000 0.00 180.00 0.00 0.0000
IC O1A CA O3A H3A 0.0000 0.00 0.00 0.00 0.0000
RESI ZERO -2.00
GROUP
ATOM CA CG2O6 0.69 ! O1
ATOM H3A HGP1 0.43 ! H3 /
ATOM O1A OG2D2 -0.76 ! O3--C (-)
ATOM O2A OG2D2 -0.76 ! \
ATOM O3A OG311 -0.60 ! O2
GROUP
ATOM C1B CG2O6 1.42
ATOM O1B OG2D2 -1.14
ATOM O2B OG2D2 -1.14
ATOM O3B OG2D2 -1.14
BOND CA O1A CA O2A CA O3A O3A H3A
BOND C1B O1B C1B O2B C1B O3B
IMPR CA O1A O2A O3A
! seed = O1A CA O2A
! required for out-of-plane vibrations
IMPR C1B O1B O2B O3B
IC O2A O1A *CA O3A 0.0000 0.00 180.00 0.00 0.0000
IC O1A CA O3A H3A 0.0000 0.00 0.00 0.00 0.0000
IC O2B O1B *C1B O3B 0.0000 0.00 180.00 0.00 0.0000
IC O3B O2B *C1B O1B 0.0000 0.00 180.00 0.00 0.0000 !redundant definition needed to enable seeding.
The FEP part in the input (forward-off.namd In-aqua) looks like this:
source ../tools/fep.tcl
alch on
alchType FEP
alchFile solvate.fep
alchCol B
alchOutFile forward-off.fepout
alchOutFreq 10
alchVdwLambdaEnd 1.0
alchElecLambdaStart 0.5
alchVdWShiftCoeff 4.0
alchDecouple off
alchEquilSteps 100
set numSteps 500
runFEP 0.0 1.0 0.05 $numSteps
Increasing the numbers of steps for equilibration and forward/backward FEP calculation didn't change the results.
While I was checking the NAMD mailing list for answers, I also noticed that people use residue patches (PRES) for alchemical perturbations with change in charge (protonation/deprotonation and pKa calculations). I didn't try this, but it makes me wonder if that would be another possibility for such calculations, where we have charged species.
I will appreciate any comments and suggestions. Please let me know if you need more information.
Thank you!
Hristina
This archive was generated by hypermail 2.1.6 : Thu Dec 31 2020 - 23:17:10 CST