Re: PMF calculation input files

From: Jérôme Hénin (jerome.henin_at_ibpc.fr)
Date: Wed Aug 24 2016 - 15:22:47 CDT

Thanks Josh for answering. I'll just add that the "colvars" command in VMD
has been renamed "cv" for consistency with NAMD, so "colvars" in VMD is
deprecated.

Jerome

On 23 August 2016 at 18:02, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov> wrote:

> The error is probably since the wallconstants you added in weren't
> removed as well? That would be my guess, since it doesn't make sense to
> have wallconstants without boundaries.
>
> Here is an example of a tcl script I wrote to calculate a collective
> variable in VMD for the last frame. More documentation is in the
> userguide (http://www.ks.uiuc.edu/Research/vmd/current/ug/node211.html)
>
> mol load psf input/system.psf namdbin input/NPT010.restart.coor
> package require pbctools
> pbc set [list 70 70 140 90 90 90]
> colvars molid 0 ; #This sets up the colvars module
> colvars configfile colvar.conf ; #This reads the configuration file, to
> see what collective variables are defined
> colvars update ; # Recalculates the values for the colvars based on the
> current frame
> set val [colvars colvar contactcount value] ; #Get the value for the
> colvar.
>
> -Josh
>
>
> On 08/23/2016 09:01 AM, Olya Kravchenko wrote:
> > Thank you!! I think Instead of defining hard boundaries for every 2A
> > window I just have to define the center of each window in each colvar
> > file (i.e. centers 135, centers 133, centers 131...etc.). Would that
> > be correct?
> >
> > Does it mean that lowerboundary and upperboundary can be commented out
> > though? When I try to comment it out there is an error that says:
> >
> > FATAl ERROR: Error in the collective variables module: exiting.
> >
> >>> For instance, when you measure the colvar in vmd as written, is the
> value really near 132?
> > How can I measure colvar in VMD?
> >
> >
> > On Mon, Aug 22, 2016 at 5:27 PM, Vermaas, Joshua
> > <Joshua.Vermaas_at_nrel.gov> wrote:
> >> Hi Olga,
> >>
> >> Window sizes are defined by how far apart you place the window centers
> >> in umbrella sampling. The harmonic potential is supposed to define the
> >> boundary quite naturally, since eventually the particle can't overcome
> >> the harmonic potential and thus has a limited sampling range.
> >> Lowerboundary and upperboundary are useful for other PMF calculation
> >> methods like metadynamics, where there isn't always a natural range, and
> >> a range needs to be imposed to keep the sampling in regions you are
> >> interested in.
> >>
> >> Your stacktrace isn't the most informative. :( If I were in your
> >> position, I'd be very systematic in figuring out where exactly the
> >> problems are coming from. For instance, is it the colvars or the
> >> restraints causing the errors? Comment out one, see if it crashes. Then
> >> try the other one. I suspect the problem lies in the colvars half (at
> >> least that's usually where something weird happens for me). For
> >> instance, when you measure the colvar in vmd as written, is the value
> >> really near 132? I'd imagine that the distance along Z between the
> >> selected atom and the center of mass of the reference should probably
> >> range between -20 and 20 if the ion is transiting through the protein,
> >> or at least be centered on 0, since distanceZ measures a relative offset
> >> between the two components, not the absolute position of an atom. I can
> >> imagine if the colvar is off by 100 Angstroms from its center, the
> >> forces would be huge (2.5 * 10000 = big number), and the stacktrace
> >> might indicate a broken assumption somewhere that the forces can't be
> >> too big.
> >>
> >> Good luck!
> >> -Josh
> >>
> >> On 08/22/2016 02:11 PM, Olya Kravchenko wrote:
> >>> Hi Joshua, thank you for replying to my questions.
> >>>
> >>> 3) Does it mean I should just remove the lowerboundary and
> >>> upperboundary options from the input file completely? How would I
> >>> define the size of the window then?
> >>>
> >>>>> What is being constrained in z and why?
> >>> I want the ion to stay within 108 and 136 along z-axis, which is the
> >>> length of the channel plus some extra distance. I thought that for
> >>> each window I have to specify physical boundaries along 1-axis as well
> >>> (i.e. such as 134 and 136).
> >>>
> >>>>> Is there a reason the width is set to something other than 1?
> >>> There is no reason, I just used the default value in the PMF tutorial.
> >>> Thanks, I've changed it to 1.
> >>>
> >>> Every time I submit my simulation it does not proceed beyond the first
> >>> step and namd quits with this in the log file:
> >>>
> >>> TCL: Running for 10000 steps
> >>> PRESSURE: 0 -1791.33 21.6158 5.78314 21.6158 -1749.04 -9.34849
> >>> -515.868 -220.473 2696.16
> >>> GPRESSURE: 0 -1753.95 21.2858 1.71399 20.5862 -1713.11 -7.46195
> >>> -514.247 -221.622 2731.84
> >>> ETITLE: TS BOND ANGLE DIHED
> >>> IMPRP ELECT VDW BOUNDARY MISC
> >>> KINETIC TOTAL TEMP POTENTIAL
> >>> TOTAL3 TEMPAVG PRESSURE GPRESSURE
> >>> VOLUME PRESSAVG GPRESSAVG
> >>>
> >>> ENERGY: 0 3903.6981 12213.6490 33309.3229
> >>> 485.1122 -1264109.2871 115215.6139 0.0000
> >>> 172.9276 0.0000 -1098808.9632 0.0000
> >>> -1098808.9632 -1098664.3685 0.0000 -281.4050
> >>> -245.0753 3125726.8315 -281.4050 -245.0753
> >>>
> >>> OPENING EXTENDED SYSTEM TRAJECTORY FILE
> >>> [0] Stack Traceback:
> >>> [0:0] _ZN12PmeRealSpace12fill_chargesEPPdS1_RiS2_PcS3_
> P11PmeParticle+0x906
> >>> [0xb987c6]
> >>> [0:1] _ZN10ComputePme6doWorkEv+0xc70 [0x9413d0]
> >>> [0:2] _ZN11WorkDistrib18messageEnqueueWorkEP7Compute+0x193
> [0xc17553]
> >>> [0:3] _ZN7Compute10patchReadyEiii+0x50 [0x57a750]
> >>> [0:4] _ZN4Sync15releaseComputesEv+0x8e [0xc04c9e]
> >>> [0:5] _ZN4Sync10PatchReadyEv+0x70 [0xc04e40]
> >>> [0:6] _ZN9HomePatch14positionsReadyEi+0x15e1 [0xadfee1]
> >>> [0:7] _ZN9Sequencer17runComputeObjectsEii+0xe6 [0xbd7ef6]
> >>> [0:8] _ZN9Sequencer9integrateEv+0x109f [0xbda06f]
> >>> [0:9] _ZN9Sequencer9algorithmEv+0x66c [0xbd568c]
> >>> [0:10] _ZN9Sequencer9threadRunEPS_+0x7 [0xbe4937]
> >>> [0:11] CthStartThread+0x1e [0xc7103e]
> >>> [0:12] +0x438f0 [0x2b59cfcca8f0]
> >>>
> >>> I've already made sure all the files are in the same directory and
> >>> found by namd, and I only have this error when I add PMF-specific
> >>> entries (colvar and restraints) into configuration files; I have also
> >>> minimized the structure and I do not see any overlaps. Now I am trying
> >>> to see if something in my PMF input files is causing this error.
> >>>
> >>> Olga
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>> On Mon, Aug 22, 2016 at 12:23 PM, Vermaas, Joshua
> >>> <Joshua.Vermaas_at_nrel.gov> wrote:
> >>>> Hi Olga,
> >>>>
> >>>> 1.) What you propose will work. Strictly speaking, the psf and pdb
> given
> >>>> just need to have the right number of atoms and give the correct
> >>>> connectivity if you are also providing a bincoordinates and
> binvelocities.
> >>>>
> >>>> 2.) The first timestep value is arbitrary, in that it will not impact
> >>>> dynamics at all. The only thing that matters to get your PMF out at
> the
> >>>> end is how you postprocess the colvars.traj files that get generated.
> >>>>
> >>>> 3.) This is where you run into problems. For conventional US sampling,
> >>>> lowerboundary and upperboundary should not be defined. You are already
> >>>> adding a harmonic potential that is supposed to keep your particle
> >>>> constrained. By defining additional wall potentials, you no longer are
> >>>> biasing the potential in the way WHAM or other postprocessing tools
> >>>> might expect. I'm not clever enough to figure out what the end result
> >>>> will be if you have half harmonic potentials on top of the underlying
> >>>> harmonic potentials, especially since the half-harmonic potentials are
> >>>> soo much stronger in your construction.
> >>>>
> >>>> 1 (again)) See the answer to 3.
> >>>>
> >>>> 2.) The harmonic potential has the form .5*k(x-center)^2, so I think
> >>>> your interpretation is correct.
> >>>>
> >>>> 3.) Yeah, that'll work, so long as the protein does not undergo large
> >>>> conformational changes.
> >>>>
> >>>>
> >>>> In reading over the files, I have two questions: What is being
> >>>> constrained in z and why? Is there a reason the width is set to
> >>>> something other than 1? Since there is only one colvar being used, its
> >>>> probably easier to think about if the force constant doesn't need to
> be
> >>>> rescaled prior to feeding it into WHAM.
> >>>>
> >>>> -Josh
> >>>>
> >>>> On 08/22/2016 09:34 AM, Olya Kravchenko wrote:
> >>>>> Hi all,
> >>>>>
> >>>>> I am trying to run US-only PMF calculations with harmonic restraints
> >>>>> for ion permeation through a three-fold pore that is part of a
> >>>>> spherical protein (according to the PMF tutorial) and getting some
> >>>>> errors that I am trying to eliminate. I would like to make sure I
> >>>>> understand everything that goes into the input files correctly.
> >>>>>
> >>>>> The reaction coordinate is z-axis that goes more or loss through the
> >>>>> middle of the channel. To create pdb files for each window I moved
> the
> >>>>> atom 2A down the channel and equilibrated the system for 300ps, then
> >>>>> took the last snapshot, and equilibrated again, thus creating 14
> >>>>> separate simulations where the ion was at a new position.
> >>>>>
> >>>>> My PMF configuration and input files are in the end of this message.
> >>>>> The questions I have are the following:
> >>>>>
> >>>>> In the configuration file:
> >>>>>
> >>>>> 1) Do I understand it correctly, that for coordinates I have to use
> >>>>> the pdb that defines the new window, and for
> >>>>> binCoordinates/extendedSystem I use the associated restart
> >>>>> equilibration files?
> >>>>>
> >>>>> 2) Is the first timestep supposed to be 0 for each separate window?
> >>>>>
> >>>>> 3) When I ran my simulations for each window there were no
> >>>>> constraints, do I still need the constraint section in the PMF
> >>>>> configuration file? I used consref ../../input/rest6.ref that is no
> >>>>> different from my original starting pdb file (after it was minimized
> >>>>> and equilibrated and before I started moving the ion through the
> >>>>> pore). Does it have to be a pdb for each new window, or the original
> >>>>> file, or no file at all?
> >>>>>
> >>>>> In the colvar input file:
> >>>>>
> >>>>> 1) Do lowerboundary and upperboundary define window size or the
> length
> >>>>> of the ion pathway? (confused because in the manual it seems to be
> the
> >>>>> length of the pathway). So if my channel starts at 136 and ends at
> >>>>> 108, I would have lowerboundary 134.0 upperboundary 136.0 for
> >>>>> window 1, lowerboundary 132.0 upperboundary 134.0 for window 2, and
> so
> >>>>> on. Is this correct?
> >>>>>
> >>>>> 2) Is "centers" a physical center of the window? So, if the window is
> >>>>> between 134 and 136 then the centers 135?
> >>>>>
> >>>>> 3) For ref I listed all CA atoms of the three chains that define the
> >>>>> pore. Does it make sense?
> >>>>>
> >>>>>
> >>>>> I'd very much appreciate comments on any of these questions. Thanks
> in advance!!
> >>>>>
> >>>>> Olga
> >>>>>
> >>>>> P.S. The input files are below:
> >>>>>
> >>>>>
> >>>>> set num 2
> >>>>>
> >>>>> #############################################################
> >>>>> ## JOB DESCRIPTION ##
> >>>>> #############################################################
> >>>>>
> >>>>> # ABF of Fe2+ in ferritin, window 1
> >>>>>
> >>>>> #############################################################
> >>>>> ## ADJUSTABLE PARAMETERS ##
> >>>>> #############################################################
> >>>>>
> >>>>> structure final.psf
> >>>>> coordinates w2.pdb
> >>>>>
> >>>>> binCoordinates w2.restart.coor
> >>>>> extendedSystem w2.restart.xsc
> >>>>>
> >>>>> set temperature 300
> >>>>> set outputname win${num}_out
> >>>>>
> >>>>> firsttimestep 0
> >>>>>
> >>>>> #############################################################
> >>>>> ## SIMULATION PARAMETERS ##
> >>>>> #############################################################
> >>>>>
> >>>>> # Input
> >>>>> paraTypeCharmm on
> >>>>> parameters ../../toppar/par_all36_prot.prm
> >>>>> parameters ../../toppar/metals/CHARMM_
> METAL/par_all22_prot_metals.inp
> >>>>> temperature $temperature
> >>>>>
> >>>>>
> >>>>> # Force-Field Parameters
> >>>>> exclude scaled1-4
> >>>>> 1-4scaling 1.0
> >>>>> cutoff 12.
> >>>>> switching on
> >>>>> switchdist 10.
> >>>>> pairlistdist 14.0
> >>>>>
> >>>>>
> >>>>> # Integrator Parameters
> >>>>> timestep 2.0 ;# 2fs/step (only if needed to finish
> quickly)
> >>>>> rigidBonds water ;# needed for 2fs steps
> >>>>> nonbondedFreq 1
> >>>>> fullElectFrequency 2
> >>>>> stepspercycle 10
> >>>>>
> >>>>>
> >>>>> # Constant Temperature Control
> >>>>> langevin on ;
> >>>>> langevinDamping 1 ;# damping coefficient (gamma) of 5/ps
> >>>>> langevinTemp $temperature
> >>>>> langevinHydrogen no ;# don't couple langevin bath to hydrogens
> >>>>>
> >>>>>
> >>>>> wrapAll on
> >>>>>
> >>>>>
> >>>>> # PME (for full-system periodic electrostatics)
> >>>>> PME yes
> >>>>> PMEGridSpacing 1.0
> >>>>>
> >>>>> # Constant Pressure Control (variable volume)
> >>>>> if {0} {
> >>>>> useGroupPressure yes ;
> >>>>> useFlexibleCell no ;
> >>>>> useConstantArea no ;
> >>>>>
> >>>>> langevinPiston on
> >>>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
> >>>>> langevinPistonPeriod 100.0
> >>>>> langevinPistonDecay 50.0
> >>>>> langevinPistonTemp $temperature
> >>>>> }
> >>>>>
> >>>>>
> >>>>> # Output
> >>>>> outputName $outputname
> >>>>>
> >>>>> restartfreq 500 ;# 500steps = every 1ps
> >>>>> dcdfreq 250
> >>>>> xstFreq 250
> >>>>> outputEnergies 100
> >>>>> outputPressure 100
> >>>>>
> >>>>> commotion yes
> >>>>>
> >>>>>
> >>>>> #############################################################
> >>>>> ## EXTRA PARAMETERS ##
> >>>>> #############################################################
> >>>>>
> >>>>> constraints on
> >>>>> consexp 2
> >>>>> consref ../../input/rest6.ref
> >>>>> conskfile ../../input/rest6.ref
> >>>>> conskcol O
> >>>>> selectConstraints on
> >>>>> selectConstrX off
> >>>>> selectConstrY off
> >>>>> selectConstrZ on
> >>>>>
> >>>>> colvars on
> >>>>> colvarsConfig US-win${num}.in
> >>>>>
> >>>>> #margin 5
> >>>>>
> >>>>> #############################################################
> >>>>> ## EXECUTION SCRIPT ##
> >>>>> #############################################################
> >>>>>
> >>>>>
> >>>>> run 150000 ;#
> >>>>>
> >>>>> Colvar input file:
> >>>>>
> >>>>> colvarsTrajFrequency 20
> >>>>>
> >>>>> colvar {
> >>>>> name Translocation
> >>>>>
> >>>>> width 0.1
> >>>>>
> >>>>> lowerboundary 132.0
> >>>>> upperboundary 134.0
> >>>>>
> >>>>> lowerwallconstant 100.0
> >>>>> upperwallconstant 100.0
> >>>>>
> >>>>> distanceZ {
> >>>>> main {
> >>>>> atomnumbers { 291491 }
> >>>>> }
> >>>>> ref {
> >>>>> atomnumbers {
> >>>>> 5 19 30 47 63 87 104 118 139 156 173 185 196 211 221 231 250 264 288
> >>>>> 305 324 338 357 372 391 412 422 433 454 470 491 510 521 538 549 570
> >>>>> 591 611 623 647 659 671 687 697 716 738 752 772 782 804 825 845 864
> >>>>> 881 898 909 926 941 956 980 995 1012 1022 1037 1059 1078 1095 1117
> >>>>> 1136 1153 1167 1184 1208 1215 1222 1246 1265 1285 1304 1321 1333 1352
> >>>>> 1369 1393 1405 1417 1428 1440 1452 1476 1491 1502 1509 1528 1542 1552
> >>>>> 1569 1584 1595 1605 1624 1641 1660 1675 1697 1711 1727 1741 1758 1769
> >>>>> 1788 1807 1822 1841 1858 1880 1899 1909 1923 1935 1957 1971 1985 1997
> >>>>> 2014 2033 2044 2056 2076 2095 2110 2124 2141 2162 2181 2195 2210 2227
> >>>>> 2243 2265 2275 2294 2316 2331 2350 2357 2369 2386 2402 2416 2430 2449
> >>>>> 2473 2495 2512 2519 2531 2543 2558 2569 2576 2595 2605 2620 2641 2660
> >>>>> 2680 2692 2714 2731 2745 2764 2774 2788 2799 2816 2832 2856 2873 2887
> >>>>> 2908 2925 2942 2954 2965 2980 2990 3000 3019 3033 3057 3074 3093 3107
> >>>>> 3126 3141 3160 3181 3191 3202 3223 3239 3260 3279 3290 3307 3318 3339
> >>>>> 3360 3380 3392 3416 3428 3440 3456 3466 3485 3507 3521 3541 3551 3573
> >>>>> 3594 3614 3633 3650 3667 3678 3695 3710 3725 3749 3764 3781 3791 3806
> >>>>> 3828 3847 3864 3886 3905 3922 3936 3953 3977 3984 3991 4015 4034 4054
> >>>>> 4073 4090 4102 4121 4138 4162 4174 4186 4197 4209 4221 4245 4260 4271
> >>>>> 4278 4297 4311 4321 4338 4353 4364 4374 4393 4410 4429 4444 4466 4480
> >>>>> 4496 4510 4527 4538 4557 4576 4591 4610 4627 4649 4668 4678 4692 4704
> >>>>> 4726 4740 4754 4766 4783 4802 4813 4825 4845 4864 4879 4893 4910 4931
> >>>>> 4950 4964 4979 4996 5012 5034 5044 5063 5085 5100 5119 5126 5138 5155
> >>>>> 5171 5185 5199 5218 5242 5264 5281 5288 5300 5312 5327 5338 5345 5364
> >>>>> 5374 5389 5410 5429 5449 5461 5483 5500 5514 5533 5543 5557 5568 5585
> >>>>> 5601 5625 5642 5656 5677 5694 5711 5723 5734 5749 5759 5769 5788 5802
> >>>>> 5826 5843 5862 5876 5895 5910 5929 5950 5960 5971 5992 6008 6029 6048
> >>>>> 6059 6076 6087 6108 6129 6149 6161 6185 6197 6209 6225 6235 6254 6276
> >>>>> 6290 6310 6320 6342 6363 6383 6402 6419 6436 6447 6464 6479 6494 6518
> >>>>> 6533 6550 6560 6575 6597 6616 6633 6655 6674 6691 6705 6722 6746 6753
> >>>>> 6760 6784 6803 6823 6842 6859 6871 6890 6907 6931 6943 6955 6966 6978
> >>>>> 6990 7014 7029 7040 7047 7066 7080 7090 7107 7122 7133 7143 7162 7179
> >>>>> 7198 7213 7235 7249 7265 7279 7296 7307 7326 7345 7360 7379 7396 7418
> >>>>> 7437 7447 7461 7473 7495 7509 7523 7535 7552 7571 7582 7594 7614 7633
> >>>>> 7648 7662 7679 7700 7719 7733 7748 7765 7781 7803 7813 7832 7854 7869
> >>>>> 7888 7895 7907 7924 7940 7954 7968 7987 8011 8033 8050 8057 8069 8081
> >>>>> 8096 8107 8114 8133 8143 8158 8179 8198 8218 8230 8252 8269 8283 8302
> >>>>> }
> >>>>> }
> >>>>> axis ( 0.0, 0.0, 1.0 )
> >>>>> }
> >>>>> }
> >>>>>
> >>>>> harmonic {
> >>>>> name Z
> >>>>> colvars Translocation
> >>>>> centers 135
> >>>>> forceConstant 0.025;# 2.5 * 0.1^2 = 0.025
> >>>>> }
> >>>>>
> >>>>>
> >>>>>
> >
>
>
>

This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:25 CST