From: Francesco Pietra (chiendarret_at_gmail.com)
Date: Thu Dec 17 2009 - 02:34:23 CST

In fact, on a rough graphic inspection (with Chimera, I have to learn
how to do that with VMD) the solvated cell dimensions are

x 147
y 176
z 138

My silly mistake is confirmed. Could you please correct the procedure
I used to measure the cell?

thanks
francesco

---------- Forwarded message ----------
From: Francesco Pietra <chiendarret_at_gmail.com>
Date: Thu, Dec 17, 2009 at 9:05 AM
Subject: Re: Fwd: vmd-l: Re: namd-l: Fwd: conf file for coarse grained
simulation
To: Peter Freddolino <petefred_at_ks.uiuc.edu>

Hello Peter:

The cell size and center came from:

--launch VMD

--Load solvated .psf and pdb

--Tk Console

--Set all [atomselect top all]

--measure center $all
  35.1028..  46.859...  28.922...

--measure minmax $all
 {-39.891...   -43.641...   -40.103...}
 {109.567...  136.141...   98.125}

(these max values are reasonable for the non-solvated protein+bilayer;
see at the end my conclusion)

--molinfo top set a 109

--molinfo top set b 136

--molinfo top set c 98

--Graphic repres.. Periodic .. select all x -x y -y z -z
to get the posted image.

In minimize conf:
# Periodic Boundary Conditions
# NOTE: Do not set the periodic cell basis if you have also
# specified an .xsc restart file!
if {1} {
cellBasisVector1    109.    0.   0.
cellBasisVector2     0.   136.   0.
cellBasisVector3     0.    0.  98.
cellOrigin         35.102874755859375 46.8595085144043 28.922983169555664
}
wrapWater           on
wrapAll             on

 #PME (for full-system periodic electrostatics)
if {1} {
PME                 yes
PMEGridSizeX       109
PMEGridSizeY       136
PMEGridSizeZ       98
}

=====
I did not take

 {-39.891...   -43.641...   -40.103...}

into account when using the Solvate plugin, where:

Boundary 5.0

Use molecule dimensions

Box padding 25 in all six boxes

Use nonstandard solvent
 cgwat.pdb
 cgwat.psf
 rbcg-2007.top
 side length 40
 name H2O
===

To create the .fix file:

--VMD ... Tk console

--mol load psf ...psf  pdb  ...pdb

--set all [atomselect top "all"]

$all set beta 0

--set fixed [atomselect top "protein or (chain O and name CHO PHO ES1
ES2 ME1 ME2 ME3 MT1 ME4 ME5 ME6 MT2)"]

$fixed set beta 1

$all writepdb ...fix

whereby POPC got 1.00 in beta column while the cg protein not (how
should it - and cg water - be named?). I set 1.00 in beta column for
the protein with a script.

>From your observations it seems to me now -as alluded to above - that
the method I used to pick up the cell dimensions did not take into
account the solvation water. If so, it was a most silly overlooking
because it was clear to me that the protein-bilayer is sized - more or
less -  109 136 98 A.

thanks
francesco

On Thu, Dec 17, 2009 at 1:49 AM, Peter Freddolino <petefred_at_ks.uiuc.edu> wrote:
> Hi Francesco,
> it is pretty apparent from the image that you posted that your periodic
> cells overlap with each other, meaning lots and lots of clashes are
> occuring. Have a look at every point where there's an interface between
> images -- you get double density of water. How did you pick your
> periodic cell size?
> Peter
>
> Francesco Pietra wrote:
>> Forgetting about the "steepest descent" procedure, probably
>> incorrectly devised, I have carefully repeated the procedure in a
>> water box.
>>
>> First I have rebuilt the protein+bilayer model with 4.0A boundary
>> between the protein and the bilayer (in previous attempts it was
>> 2.4A). The protein+bilayer could be successfully minimized
>> non-periodic, no PME, no restraint on any bead:
>> LINE MINIMIZER BRACKET: DX 9.90228e-07 7.94702e-06 DU -0.00210952
>> 0.136292 DUDX -4274.13 13.24 34182.2
>> LINE MINIMIZER REDUCING GRADIENT FROM 54035.7 TO 54.0357
>> WRITING COORDINATES TO DCD FILE AT STEP 999
>> TIMING: 1000  CPU: 36.1463, 0.0465229/step  Wall: 81.7678,
>> 0.0703083/step, 0 hours remaining, 6.966660 MB of memory in use.
>> PRESSURE: 1000 0 0 0 0 0 0 0 0 0
>> GPRESSURE: 1000 0 0 0 0 0 0 0 0 0
>> ETITLE:      TS           BOND          ANGLE          DIHED
>> IMPRP               ELECT            VDW       BOUNDARY           MISC
>>        KINETIC               TOTAL           TEMP      POTENTIAL
>>   TOTAL3        TEMPAVG
>>
>> ENERGY:    1000      1266.4280      3301.3107      1238.8463
>> 0.0000          -5076.5324    -18271.7678         0.0000
>> 0.0000         0.0000         -17541.7153         0.0000
>> -17541.7153    -17541.7153         0.0000
>>
>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 1000
>> WRITING COORDINATES TO DCD FILE AT STEP 1000
>> WRITING COORDINATES TO RESTART FILE AT STEP 1000
>> FINISHED WRITING RESTART COORDINATES
>> WRITING VELOCITIES TO RESTART FILE AT STEP 1000
>> FINISHED WRITING RESTART VELOCITIES
>> REINITIALIZING VELOCITIES AT STEP 1000 TO 300 KELVIN.
>> TCL: Setting parameter numsteps to 10
>> WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP 10
>> CLOSING EXTENDED SYSTEM TRAJECTORY FILE
>> WRITING COORDINATES TO OUTPUT FILE AT STEP 10
>> CLOSING COORDINATE DCD FILE
>> WRITING VELOCITIES TO OUTPUT FILE AT STEP 10
>> ====================================================
>>
>> WallClock: 83.922256  CPUTime: 37.538345  Memory: 6.827713 MB
>> Charmrun> Graceful exit.
>>
>> cg Water solvation of the minimized system was then carried out with
>> Boundary 5.0A and Padding 25A (15+10). In previous attempts it was 4.0
>> and 19 (15+4).
>>
>> Minimization was carried out under the same conditions as in previous
>> runs (periodic, PME, restraint on the protein and the lipid, which
>> means that no water was restrained, be that of external solvation or
>> associated with POPC)
>>
>> The simulation halted at step 542 (out of stated 1000) with (as
>> before) no warning signal. Only one processor was still active 100%
>> (top monitoring). The simulation started removing clashes (between
>> solvation waters, or the gap between cell is still not large enough?
>> please see the attached image of the periodic cells):
>>
>> OPENING EXTENDED SYSTEM TRAJECTORY FILE
>> MINIMIZER SLOWLY MOVING 18573 ATOMS WITH BAD CONTACTS DOWNHILL
>> PRESSURE: 1 -1.85929e+15 -1.99318e+15 -1.13841e+15 7.73135e+13
>> 2.58478e+15 -1.60992e+15 2.13888e+14 4.95167e+14 -7.73779e+13
>> GPRESSURE: 1 -1.85929e+15 -1.99318e+15 -1.13841e+15 7.73135e+13
>> 2.58478e+15 -1.60992e+15 2.13888e+14 4.95167e+14 -7.73779e+13
>> ENERGY:       1      1266.4376      3301.3237      1238.8804
>> 0.0000          -4957.8167 9999999999.9999         0.0000
>> 0.0000         0.0000      9999999999.9999         0.0000
>> 9999999999.9999 9999999999.9999         0.0000      9999999999.9999
>> 9999999999.9999   1452752.0000 9999999999.9999 9999999999.9999
>>
>> MINIMIZER SLOWLY MOVING 16749 ATOMS WITH BAD CONTACTS DOWNHILL
>> OPENING COORDINATE DCD FILE
>> WRITING COORDINATES TO DCD FILE AT STEP 1
>> PRESSURE: 2 -9.61887e+13 -9.14159e+13 -5.43871e+13 2.85652e+13
>> 1.3131e+14 -6.18717e+13 1.15444e+13 2.7896e+13 -6.19013e+12
>> GPRESSURE: 2 -9.61887e+13 -9.14159e+13 -5.43871e+13 2.85652e+13
>> 1.3131e+14 -6.18717e+13 1.15444e+13 2.7896e+13 -6.19013e+12
>> ENERGY:       2      1266.4376      3301.3237      1238.8804
>> 0.0000          -4957.8167 9999999999.9999         0.0000
>> 0.0000         0.0000      9999999999.9999         0.0000
>> 9999999999.9999 9999999999.9999         0.0000      9999999999.9999
>> 9999999999.9999   1452752.0000 9999999999.9999 9999999999.9999
>>
>> and continues so until step 26, than initiates conjugated gradient:
>>
>> MINIMIZER SLOWLY MOVING 1 ATOMS WITH BAD CONTACTS DOWNHILL
>> WRITING COORDINATES TO DCD FILE AT STEP 25
>> PRESSURE: 26 1.49959e+08 -7.27695e+06 -4.73455e+06 1.66896e+06
>> 1.50542e+08 -215263 -1.29055e+06 -3.77079e+06 1.51917e+08
>> GPRESSURE: 26 1.49959e+08 -7.27695e+06 -4.73455e+06 1.66896e+06
>> 1.50542e+08 -215263 -1.29055e+06 -3.77079e+06 1.51917e+08
>> ENERGY:      26      1266.4376      3301.3237      1238.8804
>> 0.0000          -4957.8167 737159432.2607         0.0000
>> 0.0000         0.0000      737160281.0857         0.0000
>> 737160281.0857 737160281.0857         0.0000      150806192.8320
>> 150806192.8320   1452752.0000 150806192.8320 150806192.8320
>>
>> MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM
>> LINE MINIMIZER REDUCING GRADIENT FROM 5.31294e+14 TO 5.31294e+11
>> WRITING COORDINATES TO DCD FILE AT STEP 26
>> PRESSURE: 27 9.27468e+07 75657.6 -1.68847e+06 1.54594e+06 9.48161e+07
>> 768612 -1.43932e+06 -1.05566e+06 9.49739e+07
>> GPRESSURE: 27 9.27468e+07 75657.6 -1.68847e+06 1.54594e+06 9.48161e+07
>> 768612 -1.43932e+06 -1.05566e+06 9.49739e+07
>> ENERGY:      27      1266.4376      3301.3237      1238.8804
>> 0.0000          -4957.8167 463067070.1446         0.0000
>> 0.0000         0.0000      463067918.9696         0.0000
>> 463067918.9696 463067918.9696         0.0000       94178916.2476
>> 94178916.2476   1452752.0000  94178916.2476  94178916.2476
>>
>> VDW energy did not improve much, until the procedure halted:
>>
>> WRITING COORDINATES TO DCD FILE AT STEP 542
>> PRESSURE: 543 8.85604e+07 1.45427e+06 -1.39861e+06 1.69754e+06
>> 9.11258e+07 841076 -1.52223e+06 -717614 9.07388e+07
>> GPRESSURE: 543 8.85604e+07 1.45427e+06 -1.39861e+06 1.69754e+06
>> 9.11258e+07 841076 -1.52223e+06 -717614 9.07388e+07
>> ENERGY:     543      1266.4376      3301.3237      1238.8804
>> 0.0000          -4957.8167 445412669.1165         0.0000
>> 0.0000         0.0000      445413517.9415         0.0000
>> 445413517.9415 445413517.9415         0.0000       90141654.9142
>> 90141654.9142   1452752.0000  90141654.9142  90141654.9142
>>
>>
>> Well, do you think that the Padding should be increased further, or
>> that I am doing some silly mistake that has nothing to do with the gap
>> between cells? I was aware of the Solvate plugin, actually I followed
>> that. I have reexamined everything without being able to detect were I
>> am wrong. Do you think worth while trying to restrain the water
>> associated with POPC?
>>
>> Thanks
>> francesco
>>
>>
>> ---------- Forwarded message ----------
>> From: Francesco Pietra <chiendarret_at_gmail.com>
>> Date: Wed, Dec 16, 2009 at 10:01 AM
>> Subject: Fwd: vmd-l: Re: namd-l: Fwd: conf file for coarse grained simulation
>> To: Axel Kohlmeyer <akohlmey_at_gmail.com>, Peter Freddolino <petefred_at_ks.uiuc.edu>
>> Cc: NAMD <namd-l_at_ks.uiuc.edu>
>>
>>
>> Attempts at minimizing under "steepest descent" ended in a crash. I
>> wonder whether a suggestion may arise from the following.
>>
>> Modifications to the conf file:
>>
>> set temperature    0
>>
>> temperature         $temperature
>>
>> langevin            off
>>
>> # Minimization
>> if {1} {
>> velocityQuenching on
>> maximumMove 1.5
>> minimize            1000
>> reinitvels          $temperature
>> }
>>
>> numsteps 10
>>
>> The log said:
>> .............
>> Info: Startup phase 7 took 0.0559781 s, 11.754 MB of memory in use
>> Info: Startup phase 8 took 0.000222921 s, 11.8803 MB of memory in use
>> Info: Finished startup at 15 s, 11.8803 MB of memory in use
>>
>> TCL: Minimizing for 1000 steps
>> ------------- Processor 2 Exiting: Caught Signal ------------
>> Signal: segmentation violation
>> Suggestion: Try running with '++debug', or linking with '-memory paranoid'.
>> [2] Stack Traceback:
>>  [0] /lib/libc.so.6 [0x7f4c538daf60]
>>  [1] _Z24sortEntries_mergeSort_v2RP12__sort_entryS1_i+0xba  [0x5ce526]
>>  [2] _ZN20ComputeNonbondedUtil32calc_pair_energy_merge_fullelectEP9nonbonded+0x351a
>>  [0x592c2a]
>>  [3] _ZN20ComputeNonbondedPair7doForceEPP8CompAtomPP11CompAtomExtPP7Results+0xae4
>>  [0x57b04e]
>>  [4] _ZN16ComputePatchPair6doWorkEv+0xa7  [0x6e21a3]
>>  [5] _ZN11WorkDistrib12enqueueWorkAEP12LocalWorkMsg+0x16  [0x92ed3c]
>>  [6] _ZN19CkIndex_WorkDistrib31_call_enqueueWorkA_LocalWorkMsgEPvP11WorkDistrib+0xf
>>  [0x92ed23]
>>  [7] CkDeliverMessageFree+0x21  [0x9c2d71]
>>  [8] _Z15_processHandlerPvP11CkCoreState+0x509  [0x9c2365]
>>  [9] CsdScheduleForever+0xa5  [0xa4bcf5]
>>  [10] CsdScheduler+0x1c  [0xa4b8f6]
>>  [11] _Z11master_initiPPc+0x280  [0x508ea0]
>>  [12] _ZN7BackEnd4initEiPPc+0x31  [0x508c19]
>>  [13] main+0x2f  [0x5045af]
>>  [14] __libc_start_main+0xe6  [0x7f4c538c71a6]
>>  [15] _ZNSt8ios_base4InitD1Ev+0x52  [0x4ff9ea]
>> Fatal error on PE 2> segmentation violation
>>
>>
>> Thanks
>> francesco
>>
>> ---------- Forwarded message ----------
>> From: Francesco Pietra <chiendarret_at_gmail.com>
>> Date: Wed, Dec 16, 2009 at 9:13 AM
>> Subject: Re: vmd-l: Re: namd-l: Fwd: conf file for coarse grained simulation
>> To: Axel Kohlmeyer <akohlmey_at_gmail.com>
>> Cc: NAMD <namd-l_at_ks.uiuc.edu>, Peter Freddolino <petefred_at_ks.uiuc.edu>
>>
>>
>> Thanks a lot for what you wrote. Solvation from scratch with Boundary
>> 4.0 and Padding 19 (15+4) did not help, or not fully. Minimization at
>> const. volume (with protein+bilayer restrained) halted at step 532,
>> out of 1000 set steps.
>>
>> The reason to posting now, before I continue to search the right
>> avenue, is that now the VDW energy has decreased conspicuously, from
>> initial
>>
>> MINIMIZER SLOWLY MOVING ATOMS WITH BAD CONTACTS DOWNHILL
>> PRESSURE: 3 1.5157e+10 1.03589e+10 4.97618e+09 9.72122e+08 1.2445e+10
>> -5.9558e+08 -2.86152e+09 -1.00705e+10 1.38611e+09
>> GPRESSURE: 3 1.5157e+10 1.03589e+10 4.97618e+09 9.72122e+08 1.2445e+10
>> -5.9558e+08 -2.86152e+09 -1.00705e+10 1.38611e+09
>> ENERGY:       3      1578.4243      3635.7565      1250.6257
>> 0.0000          -5656.7864 9999999999.9999         0.0000
>> 0.0000         0.0000      9999999999.9999         0.0000
>> 9999999999.9999 9999999999.9999         0.0000      9662690638.0470
>> 9662690638.0470   1332198.0000 9662690638.0470 9662690638.0470
>>
>> MINIMIZER SLOWLY MOVING ATOMS WITH BAD CONTACTS DOWNHILL
>> PRESSURE: 4 4.34298e+09 6.51726e+09 4.72519e+09 1.90218e+09
>> 5.44552e+09 3.08174e+09 -2.88084e+07 -4.85553e+08 9.49706e+08
>> GPRESSURE: 4 4.34298e+09 6.51726e+09 4.72519e+09 1.90218e+09
>> 5.44552e+09 3.08174e+09 -2.88084e+07 -4.85553e+08 9.49706e+08
>> ENERGY:       4      1583.3489      3638.1393      1260.5904
>> 0.0000          -5655.4262 7235356386.7152         0.0000
>> 0.0000         0.0000      7235357213.3676         0.0000
>> 7235357213.3676 7235357213.3676         0.0000      3579402107.2727
>> 3579402107.2727   1332198.0000 3579402107.2727 3579402107.2727
>>
>> to final
>>
>> WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP 530
>> WRITING COORDINATES TO DCD FILE AT STEP 530
>> WRITING COORDINATES TO RESTART FILE AT STEP 530
>> FINISHED WRITING RESTART COORDINATES
>> WRITING VELOCITIES TO RESTART FILE AT STEP 530
>> FINISHED WRITING RESTART VELOCITIES
>> PRESSURE: 531 1.53654e+07 3.29332e+06 -372766 543512 1.72077e+07
>> 1.39373e+06 1.05476e+06 1.19923e+06 1.78077e+07
>> GPRESSURE: 531 1.53654e+07 3.29332e+06 -372766 543512 1.72077e+07
>> 1.39373e+06 1.05476e+06 1.19923e+06 1.78077e+07
>> ENERGY:     531      1600.1030      3648.7438      1294.4447
>> 0.0000          -5651.9148  73889513.1000         0.0000
>> 0.0000         0.0000       73890404.4767         0.0000
>> 73890404.4767  73890404.4767         0.0000       16793614.9574
>> 16793614.9574   1332198.0000  16793614.9574  16793614.9574
>>
>> WRITING COORDINATES TO DCD FILE AT STEP 531
>> PRESSURE: 532 1.53654e+07 3.29332e+06 -372766 543512 1.72077e+07
>> 1.39373e+06 1.05476e+06 1.19923e+06 1.78077e+07
>> GPRESSURE: 532 1.53654e+07 3.29332e+06 -372766 543512 1.72077e+07
>> 1.39373e+06 1.05476e+06 1.19923e+06 1.78077e+07
>> ENERGY:     532      1600.1030      3648.7438      1294.4447
>> 0.0000          -5651.9148  73889513.1000         0.0000
>> 0.0000         0.0000       73890404.4767         0.0000
>> 73890404.4767  73890404.4767         0.0000       16793614.9574
>> 16793614.9574   1332198.0000  16793614.9574  16793614.9574
>>
>>
>> Why the minimizer found it impossible to continue is what I am
>> wondering about. I am considering:
>>
>> ---Before running conjugate gradient, run steepest descent (which is
>> the rule in Amber). I rely on Peter's 2008 suggestion: "you can get
>> something very similar by using velocityQuenching (turn on
>> velocityQuenching and then use run X to run X steps). This method
>> removes all velocity from all atoms at each step, which gives you a
>> similar effect. See
>> http://www.ks.uiuc.edu/Research/namd/2.6/ug/node29.html#8242 for more
>> details."
>>
>> ---Increase Padding or Boundary, or both.
>>
>> ---Restrain also the water pertaining to the bilayer, in order to
>> relax the external solvation water only (may be by setting 1.00 on col
>> B for POPC water).
>>
>> ---Re-building protein+bilayer with a larger boundary (present model
>> was built with 2.5A boundary between the protein and the bilayer, and
>> could be minimized with namd under non-periodic conditions. This
>> notwithstanding,  may be that a larger boundary, 4.0A or so, is
>> needed).
>>
>> That's all i can think about now.
>>
>> francesco
>>
>> On Tue, Dec 15, 2009 at 6:16 PM, Axel Kohlmeyer <akohlmey_at_gmail.com> wrote:
>>
>>> On Tue, 2009-12-15 at 17:01 +0100, Francesco Pietra wrote:
>>>
>>>> I forgot to ask: could you please suggest how roughly modify the
>>>> parameters for cg solvation? I must have misinterpreted the analysis
>>>> of inter-cell gap. In particular, the relationship between "Boundary"
>>>> and "Box padding" is not clear to me.
>>>>
>>> please have a look at the online help of the solvate command. the
>>> html file is a bit terse in that respect and should be updated.
>>> also the URL pointing to the namd tutorial is off by one node...
>>>
>>> the boundary value is the distance between the solvent and solute.
>>> the default value of 2.4 is fairly generous, but due to the increased
>>> size of water (Martini rolls 4 waters into one site) and protein
>>> side chain "atoms", stepping this up to, say, 4.0 might be a safe
>>> choice. this can be easily rationalized from applying common sense:
>>> just compare the values in an all-atom .par file to the coarse grain
>>> .par file. the r2min/2 value in the CG .par file is 2.35 whereas the
>>> corresponding AA values are between 1.3 and 1.8 with a few around 2.0.
>>>
>>> since your system was minimizing fine w/o periodicity, the default
>>> might already be mostly ok for you. these values are empirical anyways.
>>>
>>> the padding value is how much solvent outside of min/max dimensions
>>> of the solute you want to add. so, dimensions of solute _plus_
>>> padding dimensions will be the new min/max of your system. AFAIK,
>>> this does not include a "safety" (i.e. the equivalent of boundary for
>>> inter cell distance), so i would just add that value or more to your
>>> box dimensions.
>>>
>>> when visually checking for overlaps with PBC, you have to increase
>>> the diameter of your vdw spheres (i just double them for our CMM
>>> cg model).
>>>
>>> HTH,
>>>   axel.
>>>
>>>
>>>> thanks
>>>> francesco
>>>>
>>>>
>>> --
>>> Dr. Axel Kohlmeyer  akohlmey_at_gmail.com
>>> Institute for Computational Molecular Science
>>> College of Science and Technology
>>> Temple University, Philadelphia PA, USA.
>>>
>>>
>>>
>