From: Vermaas, Josh (vermaasj_at_msu.edu)
Date: Thu Mar 25 2021 - 19:19:34 CDT
Hi Carlo,
Segfaults are progress! Now comes another question. What does a force constant of 50000 mean? Your intuition is right, that your system is likely pretty well behaved on its own since it didn't crash when you ran it when no bias was effectively applied. I suspect once you math out what you are doing to your poor simulation (like what happens if there are 326 waters in the column instead of 327) that you'll find that your force constant is *very* high. My intuition is that you can probably get away with a force constant in the tens of kcal/mol/cv^2.
-Josh
On 3/25/21 4:40 PM, Carlo Guardiani wrote:
Dear Josh and Giacomo,
first of all, many thanks for your prompt reply.
Following your advice I did the following tests,
but unfortunately none of them had success.
1) Following Josh advice, I used the beta column
to mark water oxygens witha label 1 while all
other atoms were assigned a flag 0.
mol load pdb TMD_4HKR_to_6BBF_frame_0.pdb
set all [atomselect top "all"]
$all set occupancy 0
$all set beta 0
set wat [atomselect top "water and (name O)"]
$wat set occupancy 1
$wat set beta 1
$all writepdb mark_water.pdb
Using this input the simulation started running
with reasonable values of the collective variable
nwaters, but as soon as the minimization finished
the program aborted with a "Segmentation fault"
error message.
2) Since Josh told me that traditionally the beta
column had been used to store the charge, I also
tried this approach.
mol load parm7 6BBF_Ext_assembly_xleap.prmtop pdb TMD_4HKR_to_6BBF_frame_0.pdb
set natm 279195
set indmax [expr $natm - 1]
set all [atomselect top "all"]
$all set occupancy 0
$all set beta 0
for { set ind 0 } { $ind <= $indmax } { incr ind } {
set atm [atomselect top "index $ind"]
set charge [$atm get charge]
$atm set beta $charge
#$atm set beta [expr abs($charge)]
$atm delete
puts "$ind"
}
set wat [atomselect top "water and (name O)"]
$wat set occupancy 1
$all writepdb mark_water.pdb
As soon as the minimization finished I had a
series of error messages of the kind
ERROR: Constraint failure in RATTLE algorithm for atom 86370!
ERROR: Constraint failure; simulation has become unstable.
I tried using the standard fix for this kind of problems:
I increased the number of minimization steps from 100 to 1000
and I decreased the time-step from 2 fs to 1 fs. However,
also in this case after the minimization the program crashed
with a "segmentation fault" message and I also noted several
"nan" in the energy values. This behaviour could suggest a bad
starting structure. However, I remind you that when I ran the
simulations with all zeros in the beta column (which likely
made the collective variable inactive) I could run for 5 ns
without any problem, which suggests that the structure is not
so bad.
Also, I noted that when using the atom charge in the beta column
the values of the collective variable nwaters are negative (but
in absolute value correct) which does not make any sense since it
is a number of water molecules. I suppose this behaviour is due to
the fact that the charge of water oxygens is negative. I therefore
ran another test writing in the beta column the absolute value of
the charge
$atm set beta [expr abs($charge)]
This test also ended with a segmentation fault immediately after
minimization but as I expected the values of the colvar returned to
be positive and correct.
3) As a final test, since Josh told me that some people use the
beta column to store the mass value I also tried this approach.
Also in this case the job ended with a segmentation fault and lots
of "nan" in the energy values. Additionally, in this case the values
of the colvar were badly over-estimated being of the order of several
thousands instead of a few hundreds.
I really don't know what else I could try. Could you please tell me
what I need to put in the beta column to constrain the number of water
molecules in the pore of my ion channels ?
Many thanks for the help and best wishes,
Carlo
Il giorno gio 25 mar 2021 alle ore 15:52 Giacomo Fiorin <giacomo.fiorin_at_gmail.com<mailto:giacomo.fiorin_at_gmail.com>> ha scritto:
Hi Carlo, Josh is correct. A zero beta column is the most likely explanation: if not, please let us know.
When implementing the NAMD backend of the mapTotal variable (based on GridForces), I tried to keep the new functionality as consistent as possible with the use cases that Josh mentioned (custom electric fields, MDFF, etc). But that also means inheriting technicalities like this one.
FYI there is a tutorial available at https://urldefense.com/v3/__https://colvars.github.io/multi-map/multi-map.html__;!!DZ3fjg!vpzN6tjmxs60J62CsW_mkYJR6NM-jPeG0J85DFZ9RllTAh6g5RfTjhsxZ9wR6_Qbpw$ <https://urldefense.com/v3/__https://colvars.github.io/multi-map/multi-map.html__;!!HXCxUKc!kThZrhTJtxh8hNu5hx7hVLUBdtogRPwTURHse1stGPl6RRxnhAXFHj48URVZKs8$>. The NAMD manual doesn't have it linked, because I wrote it after NAMD 2.14 was released. It's geared toward variables with more than one map, but for a single map the syntax is otherwise almost identical, as I'm sure you noticed.
You may also find the latest VMD alpha build to visualize/test the variable. (Note however that VMD doesn't have GridForces, so the way map variables are defined there is slightly different).
Giacomo
On Thu, Mar 25, 2021 at 10:35 AM Vermaas, Josh <vermaasj_at_msu.edu<mailto:vermaasj_at_msu.edu>> wrote:
Hi Carlo,
What is in the beta column? You are using the occupancy column to determine if the atoms should be coupled to the potential or not. But the "weight" for each atom also matters. For traditional grid forces, when it was being used to apply an electric field, folks would copy the atomic charges to the beta column so that a carbon and hydrogen atom wouldn't feel the same attractive force. MDFF folks will use the mass of the element. The point being, is that your beta column can't be zero, otherwise it doesn't contribute to the gridforce value colvars is depending on. If you just want to count waters, I *think* the correct thing is to have the "noh and water" have a 1 in the beta column. You should then see non-zero values for the colvar.
-Josh
On 3/25/21 6:22 AM, Carlo Guardiani wrote:
Dear NAMD experts,
I am trying to run simulations with volumetric map-based
variables. At the moment I am trying to bias the number
of water molecules inside the pore of an ion channel
following the instructions in the Reference manual for
the Collective Variables Module. I built the pore occupancy
map using VMD on a preliminary trajectory
set sel [atomselect top "water and (z > -40) and (z < 25) and (sqrt(x*x+y*y) < 12)"]
volmap occupancy $sel -allframes -combine max -o Occupancy.dx
package require volutil
volutil -o Occupancy_smooth.dx -pad -smooth 0.0001 Occupancy.dx
The map, visualized with VMD, looks reasonable in that all
the pore region appears to be filled.
Following the manual, I also created a file where all water
molecules are marked in the occupancy column. Is that correct ?
Maybe I had to mark only the water molecule initially inside the
pore ? The file was generated with the following tcl script
mol load pdb TMD_4HKR_to_6BBF_frame_0.pdb
set all [atomselect top "all"]
$all set occupancy 0
set wat [atomselect top "water and (name O)"]
$wat set occupancy 1
$all writepdb mark_water.pdb
I added the following commands to the main NAMD input file
mGridForce yes
mGridForcePotFile Cavity Occupancy_smooth.dx
mGridForceFile Cavity mark_water.pdb
mGridForceCol Cavity O
mGridForceChargeCol Cavity B
mGridForceScale Cavity 0.0 0.0 0.0
colvars on
colvarsConfig ./ColvarConfig.conf
Finally, I defined in file ColvarConfig.conf a collective
variable that, I expected, should count the number of water
molecules in the cavity.
colvar {
name nwaters
mapTotal {
mapName Cavity
componentCoeff 1.0
}
}
harmonic {
colvars nwaters
centers 277 # centered around initial number of water molecules
forceConstant 50000
}
As you can see, I applied an harmonic bias to keep the number
of water molecules in the cavity close to the initial value
(277 molecules). Here comes the problem: according to the
.traj file the value of collective variable nwaters is constantly
equal to zero which is clearly wrong. The system was run for 5.0 ns
with NAMD 2.14. During the simulation the number of water molecules
inside the pore made oscillations of the order of 10% with respect
to the initial value. However, I suspect that an unbiased equilibrium
simulation would have done the same, so that this cannot be used as
a diagnostic criterium that the simulation is doing what it was
intended to do.
Could you please tell me what to do ? This is the first time I run this
kind of simulations so that some trivial mistake is highly possible.
Many thanks for your help,
Carlo Guardiani
________________________________________________________
Le informazioni contenute in questo messaggio di posta elettronica sono strettamente riservate e indirizzate esclusivamente al destinatario. Si prega di non leggere, fare copia, inoltrare a terzi o conservare tale messaggio se non si è il legittimo destinatario dello stesso. Qualora tale messaggio sia stato ricevuto per errore, si prega di restituirlo al mittente e di cancellarlo permanentemente dal proprio computer.
The information contained in this e mail message is strictly confidential and intended for the use of the addressee only. If you are not the intended recipient, please do not read, copy, forward or store it on your computer. If you have received the message in error, please forward it back to the sender and delete it permanently from your computer system.
________________________________
-- Josh Vermaas Assistant Professor, MSU-DOE Plant Research Lab and Department of Biochemisty and Molecular Biology vermaasj_at_msu.edu<mailto:vermaasj_at_msu.edu> https://urldefense.com/v3/__https://prl.natsci.msu.edu/people/faculty/josh-vermaas/__;!!DZ3fjg!vpzN6tjmxs60J62CsW_mkYJR6NM-jPeG0J85DFZ9RllTAh6g5RfTjhsxZ9wd8sVdsw$ <https://urldefense.com/v3/__https://prl.natsci.msu.edu/people/faculty/josh-vermaas/__;!!DZ3fjg!vf1OksGrvmaNdGhfe05Ths_wW2XkImtUO-ap3ohAbsjbYDwiTZu5jMMDfiw4w_hyKw$> ________________________________________________________ Le informazioni contenute in questo messaggio di posta elettronica sono strettamente riservate e indirizzate esclusivamente al destinatario. Si prega di non leggere, fare copia, inoltrare a terzi o conservare tale messaggio se non si è il legittimo destinatario dello stesso. Qualora tale messaggio sia stato ricevuto per errore, si prega di restituirlo al mittente e di cancellarlo permanentemente dal proprio computer. The information contained in this e mail message is strictly confidential and intended for the use of the addressee only. If you are not the intended recipient, please do not read, copy, forward or store it on your computer. If you have received the message in error, please forward it back to the sender and delete it permanently from your computer system. ________________________________ -- Josh Vermaas Assistant Professor, MSU-DOE Plant Research Lab and Department of Biochemisty and Molecular Biology vermaasj_at_msu.edu<mailto:vermaasj_at_msu.edu> https://urldefense.com/v3/__https://prl.natsci.msu.edu/people/faculty/josh-vermaas/__;!!DZ3fjg!vpzN6tjmxs60J62CsW_mkYJR6NM-jPeG0J85DFZ9RllTAh6g5RfTjhsxZ9wd8sVdsw$
This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:11 CST