Fw: [gmx-users] position restraints
Justin Lemkul
jalemkul at vt.edu
Tue Mar 26 18:25:50 CET 2013
On Tue, Mar 26, 2013 at 1:01 PM, Shima Arasteh
<shima_arasteh2001 at yahoo.com>wrote:
>
>
>
> Have a look at processed topology file here please; I see that position
> restraints are brought after chain_A but not brought after chain_B.
>
> With these settings:
> ; Include chain topologies
> #ifdef POSRES
> #include "topol_Protein_chain_A.itp"
> #include "protein_chain_A_posre.itp"
> #endif
> #ifdef POSRES
> #include "topol_Protein_chain_B.itp"
> #include "protein_chain_B_posre.itp"
> #endif
>
>
The above approach is incorrect. The inclusion of protein topologies is
dependent upon using position restraints? That's certainly not right,
especially if you ever want to run a simulation without restraints. The
following is the correct approach:
; Include chain topologies
#include "topol_Protein_chain_A.itp"
#ifdef POSRES
#include "protein_chain_A_posre.itp"
#endif
#include "topol_Protein_chain_B.itp"
#ifdef POSRES
#include "protein_chain_B_posre.itp"
#endif
Also adding "define = -DPOSRES_LIPID -DPOSRES ", I get this processed.top:
>
>
> #grompp -f nvt.mdp -c minim.gro -p topol.top -n index.ndx -o nvt.tpr -pp
>
> THIS IS THE PROCESSED TOPOLOGY:
> ; File 'topol.top' was generated
> ; By user: shima (1000)
> ; On host: linux-cbyo.site
> ; At date: Wed Dec 12 12:17:51 2012
> ;
> ; This is a standalone topology file
> ;
> ; It was generated using program:
> ; pdb2gmx - VERSION 4.5.5
> ;
> ; Command line was:
> ; pdb2gmx -f dimer-rotated.pdb -o dimer-processed.pdb -ter -water tip3p
> ;
> ; Force field was read from current directory or a relative path - path
> added.
> ;
>
> ; Include forcefield parameters
> ; Conversion of CHARMM36 parameters to GROMACS format by Thomas Piggot,
> July 2010.
> ; Also added some parameters from later CHARMM27 versions, such as those
> for the PS headgroup.
> ;
> ; If you use these parameters please check out the forcefield.doc for
> papers to cite
>
> [ defaults ]
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 2 yes 1.0 1.0
>
> ;
> ; Added new and changed old atomtypes and pairtypes (OSL, OSLP, HBL and
> CCL) - Thomas Pigot July 2010
> ;
> [ atomtypes ]
> ;name at.num mass charge ptype sigma epsilon
> C 6 12.01100 0.51 A 0.356359487256 0.46024
> CA 6 12.01100 -0.115 A 0.355005321205 0.29288
> CC 6 12.01100 0.62 A 0.356359487256 0.29288
> CD 6 12.01100 0.000 A 0.356359487256 0.29288 ; partial
> charge def not found
> CE1 6 12.01100 0.000 A 0.372395664183 0.284512 ; partial
> charge def not found
> .
> .
> .
> .
> X NN1 CN1A X 9 180.0 10.46 2
> X NN2 CN3B X 9 180.0 4.184 2
> X CN3 CN3C X 9 180.0 0.4184 2
> X NN2 CN3C X 9 180.0 0.4184 2
> X ON4 P3 X 9 0.0 1.2552 3
>
> [ dihedraltypes ]
> ; i j k l func q0 cq
> NN2B CN4 CN5 HN2 2 0.0 58.576
> NN2G CN4 CN1 HN2 2 0.0 6.6944
> NN1 CN2 HN1 HN1 2 0.0 50.208
> CN1 NN2G CN5G ON1 2 0.0 753.12
> CN1T NN2B NN2U ON1 2 0.0 920.48
> CN1 NN2U CN3T ON1 2 0.0 753.12
> CN2 NN3G NN2G NN1 2 0.0 334.72
> CN2 NN3A CN5 NN1 2 0.0 334.72
> CN2 NN3 CN3 NN1 2 0.0 502.08
> CN4 NN2G NN3I HN3 2 0.0 326.352
> CN3 CN3C CN8 HN6 2 0.0 125.52
> HN2 CN3 CN3B NN2 2 0.0 418.4
> HN1 HN1 CN1A NN1 2 0.0 -41.84
> HN8 CN3 CN3 CN8 2 0.0 150.624
> HR1 NR1 NR2 CPH2 2 0.00 4.184
> HR1 NR2 NR1 CPH2 2 0.00 4.184
> HR3 CPH1 NR1 CPH1 2 0.00 4.184
> HR3 CPH1 NR2 CPH1 2 0.00 4.184
> HR3 NR1 CPH1 CPH1 2 0.00 4.184
> HR3 NR2 CPH1 CPH1 2 0.00 4.184
> NR1 CPH1 CPH2 CN7B 2 0.00 5.0208
> NR1 CPH2 CPH1 CN7B 2 0.00 5.0208
> HN2 X X NN2 2 0.0 8.368
> HN1 X X NN1 2 0.0 33.472
> CN1 X X ON1 2 0.0 753.12
> CN1T X X ON1 2 0.0 753.12
> CN1 X X ON1C 2 0.0 669.44
> CN2 X X NN1 2 0.0 753.12
> CN9 X X CN3T 2 0.0 117.152
> HN3B X X CN3 2 0.0 125.52
> HN3B X X CN3A 2 0.0 108.784
> HN3B X X CN3B 2 0.0 108.784
> ON1 X X CN1A 2 0.0 334.72
> HN3 X X CN3C 2 0.0 443.504
> HN6 X X CN3C 2 0.0 443.504
>
>
> ; Include chain topologies
> ;
> ; File 'topol_Protein_chain_A.itp' was generated
> ; By user: shima (1000)
> ; On host: linux-cbyo.site
> ; At date: Wed Dec 12 12:17:55 2012
> ;
> ; This is a include topology file
> ;
> ; It was generated using program:
> ; pdb2gmx - VERSION 4.5.5
> ;
> ; Command line was:
> ; pdb2gmx -f dimer-rotated.pdb -o dimer-processed.pdb -ter -water tip3p
> ;
> ; Force field was read from current directory or a relative path - path
> added.
> ;
>
> [ moleculetype ]
> ; Name nrexcl
> Protein_chain_A 3
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass
> typeB chargeB massB
> ; residue 1 FVAL rtp FVAL q 0.0
> 1 C 1 FVAL CN 1 0.357 12.011 ;
> qtot 0.357
> 2 HA 1 FVAL H1 2 0.1 1.008 ;
> qtot 0.457
> 3 O 1 FVAL ON 3 -0.51 15.999 ;
> qtot -0.053
> 4 NH1 1 FVAL N 4 -0.423 14.007 ;
> qtot -0.476
> 5 H 1 FVAL HN 5 0.333 1.008 ;
> qtot -0.143
> .
> .
> .
> 359 OC 24 GLY OT1 359 -0.67 15.999 ;
> qtot 1.67
> 360 OC 24 GLY OT2 360 -0.67 15.999 ;
> qtot 1
>
> [ bonds ]
> ; ai aj funct c0 c1 c2 c3
> 1 2 1
> 1 3 1
> 1 4 1
> .
> .
> .
> .
> 357 360 1
>
> [ angles ]
> ; ai aj ak funct c0 c1
> c2 c3
> 2 1 3 5
> 2 1 4 5
> .
> .
> .
> .
> 355 358 360 5
> 359 358 360 5
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3 c4 c5
> 2 1 4 5 9
> 2 1 4 6 9
> 3 1 4 5 9
> .
> .
> .
> 357 355 358 359 9
> 357 355 358 360 9
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3
> 1 4 3 2 2
> 20 18 22 21 2
> 29 22 31 30 2
> 31 29 33 32 2
> 48 33 50 49 2
> 50 48 52 51 2
> 55 52 57 56 2
> 57 55 59 58 2
> 74 59 76 75 2
> 76 74 78 77 2
> 85 78 87 86 2
> 87 85 89 88 2
> 104 89 106 105 2
> 106 104 108 107 2
> 114 108 116 115 2
> 116 114 118 117 2
> 134 118 136 135 2
> 136 134 138 137 2
> 145 138 147 146 2
> 147 145 149 148 2
> 161 149 163 162 2
> 163 161 165 164 2
> 171 165 173 172 2
> 173 171 175 174 2
> 187 175 189 188 2
> 189 187 191 190 2
> 198 191 200 199 2
> 200 198 202 201 2
> 217 202 219 218 2
> 219 217 221 220 2
> 227 221 229 228 2
> 229 227 231 230 2
> 251 231 253 252 2
> 253 251 255 254 2
> 262 255 264 263 2
> 264 262 266 265 2
> 282 266 284 283 2
> 284 282 286 285 2
> 292 286 294 293 2
> 294 292 296 295 2
> 307 313 310 309 2
> 316 296 318 317 2
> 318 316 320 319 2
> 327 320 329 328 2
> 329 327 331 330 2
> 342 348 345 344 2
> 351 331 353 352 2
> 353 351 355 354 2
> 358 355 360 359 2
>
> [ cmap ]
> ; ai aj ak al am funct
> 18 20 22 29 31 1
> 29 31 33 48 50 1
> 48 50 52 55 57 1
> .
> .
> .
> 292 294 296 316 318 1
> 316 318 320 327 329 1
> 327 329 331 351 353 1
>
> ; Include Position restraint file
> ; In this topology include file, you will find position restraint
> ; entries for all the heavy atoms in your original pdb file.
> ; This means that all the protons which were added by pdb2gmx are
> ; not restrained.
>
>
Here's where you start to have problems. This directive was definitely
created by pdb2gmx, which by default restrains only heavy atoms.
> [ position_restraints ]
> ; atom type fx fy fz
> 1 1 1000 1000 1000
> 3 1 1000 1000 1000
> 4 1 1000 1000 1000
> 6 1 1000 1000 1000
>
Now here is another restraint directive that clearly you have created,
which contradicts the directive above. You should only ever have one [
position_restraints ] directive in your topology. I don't know how this
would have arisen, unless you have #included a position restraint file
within the protein chain topology, in which case you're actually #including
two restraint definitions. Don't do that.
> ; position restraints for a_1-360 of Protein
>
> [ position_restraints ]
> ; i funct fcx fcy fcz
> 1 1 100000 100000 100000
> 2 1 100000 100000 100000
> 3 1 100000 100000 100000
> 4 1 100000 100000 100000
>
>
<snip>
[ moleculetype ]
> ; Name nrexcl
> Protein_chain_B 3
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass
> typeB chargeB massB
> ; residue 1 FVAL rtp FVAL q 0.0
> 1 C 1 FVAL CN 1 0.357 12.011 ;
> qtot 0.357
> 2 HA 1 FVAL H1 2 0.1 1.008 ;
> qtot 0.457
> 3 O 1 FVAL ON 3 -0.51 15.999 ;
> qtot -0.053
> 4 NH1 1 FVAL N 4 -0.423 14.007 ;
> qtot -0.476
> .
> .
> .
> .
> 359 OC 24 GLY OT1 359 -0.67 15.999 ;
> qtot 1.67
> 360 OC 24 GLY OT2 360 -0.67 15.999 ;
> qtot 1
>
> [ bonds ]
> ; ai aj funct c0 c1 c2 c3
> 1 2 1
> 1 3 1
> 1 4 1
> 4 5 1
> 4 6 1
> 6 7 1
> .
> .
> .
> .
>
>
>
> [ pairs ]
> ; ai aj funct c0 c1 c2 c3
> 1 7 1
> 1 8 1
> 1 18 1
> 2 5 1
> .
> .
> .
>
>
>
> [ angles ]
> ; ai aj ak funct c0 c1
> c2 c3
> 2 1 3 5
> 2 1 4 5
> 3 1 4 5
> 1 4 5 5
> .
> .
> .
> .
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3 c4 c5
> 2 1 4 5 9
> 2 1 4 6 9
> 3 1 4 5 9
> 3 1 4 6 9
> .
> .
> .
>
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3
> 1 4 3 2 2
> 20 18 22 21 2
> 29 22 31 30 2
> 31 29 33 32 2
> .
> .
> .
>
> [ cmap ]
> ; ai aj ak al am funct
> 18 20 22 29 31 1
> 29 31 33 48 50 1
> .
> .
> .
>
>
> ; Include Position restraint file
>
> ; position restraints for a_360-720 of Protein
>
>
Again, this is clearly wrong since you start restraining from atom 360 and
onward, so this must have come from your incorrect usage of genrestr
before. Restraints can't be applied using these atom numbers and should
have triggered a fatal error.
> [ position_restraints ]
> ; i funct fcx fcy fcz
> 360 1 100000 100000 100000
> 361 1 100000 100000 100000
>
>
My suggestion would be to revert to your original topology, do not #include
any custom restraints to your protein, because (1) you almost certainly do
not need them and (2) you have not created them correctly. The only
modification you need to make is with respect to the lipid restraints,
which are #included correctly based on your previous messages.
-Justin
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540)
231-9080http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list