CHARMM Development Project Forums User Discussion & Questions Parameter Set Discussion Generalized Born and carbohydrates

I have another problem using the Monte Carlo module in CHARMM. There must be errors while computing the energy because the Delta-E column in the energy output contains values between 0 and several thousands. Delta-E is the difference between the MC running total and the current total (according to mc.doc). I have set "INBFrq -1 IMGFrq -1", but I have also played with different values between 0 and 100. I am using GBMV (method II) for solvation with nonbond options: 16/14/12. However, I experience the same amount of errors with any other implicit solvation model (including none) or cutoffs. The errors are bigger with larger Monte Carlo moves.

The system consists of a rigid protein placed 20A away from a silicate surface. I want to find the orientation with the lowest energy. At the beginning of the Monte Carlo run, the protein spins around fast and the values in Delta-E are large, but once the protein is in contact with the surface the values in Delta-E are smaller and become 0 once the protein stops moving.

The example input file says:

! IECHeck is the frequency of comparing the MC running total of the energy
! with a call to the standard ENERGY routine. The difference is printed in
! the Delta-E column of the "MC E>" table.
...
! Small errors typically derive from atoms previously outside CUTNB.
! Large errors typically derive from atoms previously outside CUTIM.

Increasing CUTIM does not help. I am using the input script in the silicates directory to set up the images and I have set Z=200 to make sure that the protein is completely contained in the main cell. This is part of my input file:

! Manipulation of images
update inbf 0 ihbf 0 imgfrq 100 cutim 18

update inbfrq -1 imgfrq -1 ihbfrq 0 -
atom CDIE eps 1 cutnb 16 ctofnb 14 ctonnb 12 switch vswitch

! Turn on GBMV method II and add SASA parameter
GBMV EPSILON 80 BUFR 0.2 MEM 20 CUTA 20 ALFRQ 1 -
GEOM BETA -12 P1 0.45 P2 1.25 P3 0.65 P6 8.0 -
CORR 1 SHIFT -0.1 SLOPE 0.9 WTYP 1 NPHI 5 -
FAST 1 SGBFRQ 4 SXD 0.3 -
SA 0.01500 SB 0.90000

! Rigid body translations of the protein
MOVE ADD MVTP RTRN BYALl nsele 1 sele segid PROT end WEIGht 1.0 DMAX 1 LABEL WTRN

! Rigid body rotations
MOVE ADD MVTP RROT BYALl nsele 1 sele segid PROT end WEIGht 1.0 DMAX 25.0 LABEL WROT

! Link translations and rotations of water so they are done together.
! Only WTRN will appear in MCSTAT.
MOVE LINK LAB1 WTRN LAB2 WROT

MC TEMPerature 300.00 NSTEps 50000 ISEEd 518282 -
IDOMcfrq 5 INBFrq -1 IMGFrq -1 IECHeck 100 -
IUNCrd 34 NSAVc 100 IACCept 0

This is part of the output:

MC E ENR: Eval# ENERgy Delta-E GRMS
MC E INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
MC E CROSS: CMAPs
MC E EXTERN: VDWaals ELEC HBONds ASP USER
MC E IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
MC E PBEQ: PBnp PBelec GBEnr GSBP
MC E PRESS: VIRE VIRI PRESSE PRESSI VOLUme
---------- --------- --------- --------- --------- ---------
MC E> 1 -604102.27176 -2176.24541 16.04694
MC E INTERN> 10207.12559 28605.58078 239.37157 6025.11303 104.29873
MC E CROSS> -237.18413
MC E EXTERN> -34290.43971-593589.79137 0.00000 598.20879 0.00000
MC E IMAGES> -2324.31511 -38728.29427 0.00000 0.00000 0.00000
MC E PBEQ> 0.00000 0.00000 19288.05434 0.00000
MC E PRESS> 1142575.9532 -355619.3542 0.0000 0.0000 0.0000
---------- --------- --------- --------- --------- ---------

Any help would be greatly appreciated! I have played around with all possible parameters but can't figure it out. Thanks!

In an explicit water MC simulation, the water molecules should also be subjected to Monte Carlo moves. Under these circumstances, the original Metropolis Monte Carlo method is expressly designed to yield a canonical ensemble, and variants for other ensembles are available in CHARMM. Over the last 50 years, MC has extensively been used for studies of bulk liquids and solute-solvent systems.

The catch is that the magnitude of the moves must be sufficiently small to get a reasonable acceptance rate; I don't know how feasible it is to achieve your goal with such small moves. It would be like doing an explicit water MD simulation, only faster (but not necessarily fast enough for your purpose).

Edit:
Your setup is unusual to me, but wouldn't all of that that be expected behavior?

Unless you're annealing (which I don't see in your input), this is not supposed to happen. I speculate you're ending up in a local minimum and your move size is so large that every step is rejected by the Metropolis criterion. This is typically not desirable because it prevents you from sampling other minima. If so, you probably want to increase the MC temperature in order to cross some barriers; as you're keeping your protein frozen and you don't seem to be aiming to recreate any kind of physical ensemble, the MC temperature can be thought of as a parameter that determines the acceptance rate.

Trying to gain insight into the problem by reading might in this case be more effective than a manual MC search of parameter space...

Kenno,


Why is there a difference between the MC running energy and the energy calculated by calling the ENERGY routine? According to the test file "c29test/mctest03.inp" Delta-E should always be 0. Quote from the test file:

! IECHeck is the frequency of comparing the MC running total of the energy
! with a call to the standard ENERGY routine. The difference is printed
! in the Delta-E column of the "MC E>" table. Given the cutoffs and update
! frequencies selected for this example, a non-zero Delta-E value for
! an IECHEck indexed larger than zero is indicative of a coding error.

This is for example what I get for the energies from a MC run (by grepping MC E>):

Step ENERgy Delta-E
0 -604687 0.00000
1 -604020 -668.01644
2 -604370 -527.25194
3 -604685 -1839.71082
4 -604680 -54.34625
5 -604680 0.00000
6 -604680 0.00000
7 -604682 -10.38353
8 -604681 -2.04802
9 -604681 0.00000
10 -604681 0.00000

I need to add that in the energy output by the MC module, Delta-E is not the difference between the current and the previous step (like for example in a minimization). It is the difference between the MC running total of the energy and a call to the standard ENERGY routine.

You are right that reading is the first thing to do to solve a problem. However, my problem is not how to run a Monte Carlo simulation. I'm trying to understand why there is a discrepancy between the total energy computed by the MC module and the output of the standard ENERGY routine (reported in the Delta-E column).

But I like your suggestion to increase the temperature to improve sampling. Thanks!

Note that in some cases it may be necessary to contact the code developer directly for such a detailed technical question; not all developers are well represented in these forums. Then again, this thread is not in the MC forum ...

In that case, I misinterpreted the CHARMM documentation. In my defense, mc.doc could be organized better - the meaning of the term "MC running total" is elusive unless one reads the section "shortcomings":Wait a second... I don't see GBMV on this list... looks like I made the right quip for the wrong reason.

So to solve the problem, either you have to switch solvent model, or instruct MC to call the main ENERGY routine at every step. It doesn't seem to be explicitly mentioned in the documentation how to do this (don't take my word for it), unless simply setting "IECHeck 1" does the trick...

Edit: of course, there's also the option of giving up on CHARMM and trying something entirely different.

Hi Kenno,

Thanks for your help. I really appreciate. I have e-mailed Aaron Dinner about this issue. He wrote back that MC does not call GBMV, only ACE and EEF1. If I want to use GBMV I would have to write an interface similar to mcace.src. EEF1 could be an option but it also uses the RDIE approximation for the electrostatics and I'm trying to find something more accurate. I am currently looking into ACE. There is a acepar22_prot.inp which contains atom volumes but of course not for silicates. So, I would have to understand how to create atom volumes for silicates.

However, in the long run I do want to use a different code than CHARMM because I want something faster and for other reasons. Thanks for the hint about the alternative software. I will look into that.

The reason why I started with CHARMM is that I wanted to evaluate different implicit solvation models and understand whether the use of RDIE to mimic the electrostatic screening of solvent is really too rough for my problem. RDIE in combination with MC in CHARMM was used in the past to predict protein orientations on "Self Assembled Monolayers": http://jcp.aip.org/resource/1/jcpsa6/v125/i21/p214704_s1?view=fulltext

Assuming that "IECHeck 1" doesn't do the trick and that I didn't miss something in the documentation, I don't understand why they don't have an option to call the full main ENERGY routine at every step. It should be trivial to code. Yes, enabling it will slow down the simulation a lot, but also opens a vast number of possibilities that are currently inaccessible for anyone unwilling or unable to modify the source. I thought giving the user an overwhelming number of possibilities is what CHARMM is all about?!

The problem is that IECHeck applies only to those steps that are not rejected. I might try to figure out how to force mcener to call ENERGY in the code, or alternatively use ACE as an implicit solvent model.

I wonder whether there is another way through which I can use CHARMM to randomly sample orientations of a rigid protein and then manually call ENERGY to find the energy minimum. Not as elegant as an MC simulation though. Well, I guess VMD would do the trick to produce such orientations. NAMD now supports a GB based implicit solvent model, but it does not have a SASA based nonpolar term. (And now the discussion is leading outside CHARMM. )