Re: PMF calculation input files

From: Olya Kravchenko (ovkrav_at_gmail.com)
Date: Mon Aug 22 2016 - 15:10:57 CDT

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:24 CST