CHARMM Development Project
Posted By: gianluca Generalized Born and carbohydrates - 04/24/11 02:46 AM
Is it possible to simulate a complex between a protein and carbohydrates (e.g., mannose) with a Generalized Born implicit solvation model like GBMV or GBSW?

Thanks!

Gianluca
Posted By: rmv Re: Generalized Born and carbohydrates - 04/24/11 02:57 AM
Possible, but perhaps unwise without extensive testing, esp. if the sugars are exposed to solvent; the new carbohydrate FF was not available when the GB* methods were developed.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/26/11 12:28 AM
The documentation of GBSW talks about a self-consistent GBSW force field:

*** SPECIAL NOTE ***

1. A self-consistent GBSW force field is optimal for peptide and protein simulations (as of 2005).

Can this self-consistent force field be used also with sugars or is it better to use the standard force field?

Also, I have another question concerning the cutoffs. The examples use: ctonnb 16 ctofnb 16 cutnb 20
Is it normal to use these cutoffs with Generalized Born? They are larger than the defaults of PARAM22.

Thank you very much!

Gianluca
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/26/11 12:49 AM
Another question. Does GBMV take into account the non-polar contributions? Does this need to be turned on by specifying the SA parameter? If SA is 0 (default), does it mean that the non-polar contributions are not taken into account?
Posted By: rmv Re: Generalized Born and carbohydrates - 04/26/11 02:27 AM
With GB* methods, one should use the custom cutoffs indicated in gb*.doc

The special GBSW force field was designed for peptides/proteins

The citations in the gb*.doc files may be useful.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/26/11 02:57 AM
Thanks! I went through the citations, but I'm still confused about the SASA module in GBMV. For example, in Lee et al 2003 it discusses how SASA is calculated, but it is not clear whether this is also used to calculate the non-polar contributions. I tend to believe that GBMV calculates only the electrostatics.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/26/11 03:16 AM
Gianluca: what nonpolar contributions do you mean? Interactions between a solute and "solvent"? Then no, there's no such thing, as far as I know. There is only electrostatic solvation energy.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/26/11 03:24 AM
What I mean is a SASA term that approximates the nonpolar contribution to solvation. This is explained in gbsw.doc:

"The GBSW module provides the (electrostatic + nonpolar) solvation
energy and forces. A Generalized Born method is used for the
electrostatic part and the solvent-exposed surface area for the
nonpolar part with a phenomenological surface tension coefficient."

It seems that GBSW has such a term whereas GBMV does not. Solvation energy is the sum of electrostatic (approximated by the GB equation) and nonpolar (approximated by the SASA) solvation energy.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/26/11 03:51 AM
oh, sorry, I misunderstood.

GBMV also has SASA, we use constants

SA = 0.01500 kcal/(mol*A**2) and SB = 0.90000 kcal/mol.

I've looked at the gbmv.doc, I'm not sure why the defaults are 0.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/26/11 06:14 AM
Thanks. That's ok. I wonder why it is called "surface area coefficient" instead of "surface tension coefficient" like in GBSW. Also, what is SB (surface area constant)?

Do you mind posting an input file which contains the GBMV setup and dynamics, if you have it handy?

Thanks!
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/26/11 11:10 PM
Quote:
I wonder why it is called "surface area coefficient" instead of "surface tension coefficient" like in GBSW.


I'm guessing different people wrote the code and documentation? smile

Quote:
Also, what is SB (surface area constant)?


I've looked at the source, it's just an offset that's added to the SASA energy (that's reported as ASP in the log). Seems like just one more fudge factor wink

Quote:
Do you mind posting an input file which contains the GBMV setup and dynamics, if you have it handy?


sure, here's what we did for 20-residue protein temperature replica exchange:

GBMV TOL 1E-10 MEM 20 CUTA 20 DN 1.0 -
BUFR 0.2 -
EPSILON 80 -
BETA -12 -
SHIFT -0.141851851851852 -
ESHIFT 0 -
SLOPE 1 -
LAMBDA1 0.5 -
P1 0.45 -
P2 1.25 -
P3 0.65 -
P6 8 -
ONX 1.9 -
OFFX 2.1 -
CORR 2 -
A1 0.3255 -
A2 0 -
A3 1.07603305785124 -
SON 1.2 SOFF 1.5 -
FAST 1 -
SGBFRQ 4 -
SXD 0.3 -
WTYP 1 NPHI 5 SA 0.015 SB 0.9

Also, as Rick noted above, you need to use cutoff scheme and the switching function suggested in gbmv.doc (switch). You can also use larger or "no cutoff" for small molecules.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/27/11 01:02 AM
And yeah, for dynamics - we often use leap/Langevin with fbeta of 10, since it imitates the collisions with water molecules that otherwise would be absent.
Posted By: rmv Re: Generalized Born and carbohydrates - 04/27/11 01:34 AM
Note that the dynamics timescale with Langevin dynamics is somewhat arbitrary. An FBETA value more like 50 gives a motional timescale closer to that in water. If conformational sampling is the goal, an FBETA value of 1 or 2 is best, as the rate of torsional barrier crossing is optimal in that regime.

Langevin Dynamics of Peptides: the Frictional Dependence of Isomerization Rates of N- Acetylalanyl-N'-Methylamide,
Richard J. Loncharich, Bernard R. Brooks and Richard W. Pastor,
Biopolymers 32, 523-535 (1992).
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 07:46 AM
Thank you for all your suggestions and advice.

I'm having trouble using the recommended cutoffs. With "cutnb 21 ctofnb 18 ctonnb 16" I get a segmentation fault. It works only if I reduce it down to 16/14/12. I'm using the xxlarge version. I'm trying to compile the huge version but I get "relocation truncated to fit: R_X86_64_32" errors during linking. I wonder whether the segmentation fault is due to the size of my system: 12000 atoms.

I'm also wondering why in the case of GBSW ctofnb=ctonnb=16 with switch? In this case, we are simply truncating the potential.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 07:49 AM
apredeus, do you have any experience with the Nose thermostat? I see that it is used in c31test/gbsw.inp.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/27/11 03:41 PM
Quote:
I'm having trouble using the recommended cutoffs. With "cutnb 21 ctofnb 18 ctonnb 16" I get a segmentation fault. It works only if I reduce it down to 16/14/12.


sounds like nonbonded arrays are becoming too large for your system. Which version of CHARMM are you using? If you can get your hands on C36a3 and later versions, they handle memory allocation better.

Quote:
I'm using the xxlarge version. I'm trying to compile the huge version but I get "relocation truncated to fit: R_X86_64_32" errors during linking. I wonder whether the segmentation fault is due to the size of my system: 12000 atoms.


Building HUGE version of CHARMM is a bit of a pain and requires removing some modules from pref.dat (or using -mcmodel=medium, but somehow it never helped on clusters I've been working on).

BUT! The thing is, actually, you probably want a smaller compilation. The sizes mean the following:

Quote:
! HUGE version - 1,000,000 atoms
! XXLARGE version - ~360,000 atoms
! XLARGE version - ~240,000 atoms
! LARGE version - ~60,000 atoms
! MEDIUM version - ~25,000 atoms
! SMALL version - ~6,000 atoms
! XSMALL version - ~2,000 atoms


as you can see, you're easily within the limits of MEDIUM. If you compile it as MEDIUM, you'll have less space used by all the main arrays, and you'll have more space left for your non-bonded and GB arrays. 12000 atoms is big but it's not outside of limits for GBMV/GBSW, you should be able to run this system.

Quote:

I'm also wondering why in the case of GBSW ctofnb=ctonnb=16 with switch? In this case, we are simply truncating the potential.


yeah, 16 is a very generous cutoff and whatever you do doesn't really matter anyway. If you truncate it at 16 you'll have some derivative discontinuity at 16 but other than that you shouldn't get any big force errors etc. You can search these forums, I think Prof. Charles Brooks answered this question before in a bit more detail.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/27/11 04:17 PM
Quote:
apredeus, do you have any experience with the Nose thermostat? I see that it is used in c31test/gbsw.inp.


yes, we've used it a lot. Usually for simple MD simulations Langevin is preferred in case of implicit solvent, but when you do temperature replica exchange you want to have very good temperature control so that your system is well equilibrated after increase/decrease of temperature and the following e.g. 500 steps of simulation. Langevin tends to lag behind in these cases, in terms of both final and average temperatures of the system.

There is an opinion that NOSE/VVER is problematic. I did not see any test cases that would unambiguously prove it, furthermore, in whatever tests we did it behaves very similarly to a newer implementation, TPCONTROL/VV2. But TPCONTROL is a lot more transparent while NOSE/VVER is a bit of a black box.

TPCONTROL also works with GBMV, so it's your call what to use.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 05:04 PM
Quote:
sounds like nonbonded arrays are becoming too large for your system. Which version of CHARMM are you using?

I'm using version c35b1.

Quote:
as you can see, you're easily within the limits of MEDIUM. If you compile it as MEDIUM, you'll have less space used by all the main arrays, and you'll have more space left for your non-bonded and GB arrays.

I have just tried with medium. The last line printed is "VEHEAP> Expanding heap size by 10305536 words."
and then it exits with "segmentation fault".

In the future, I want to use GB also for dynamics. Right now I'm using it for Monte Carlo simulations where I'm predicting the orientation of a protein on a quartz surface. The protein and surface are rigid. Thus I wonder whether a shorter cutoff of 16/14/14 could be enough. In the past a distance dependent dielectric was used, but I want to be more accurate and use GB for solvation electrostatics and a SASA term for the nonpolar.

I want to try to get the huge version to compile. Which features from pref.dat do you recommend turning off? I tried compiling with the LITE keyword. It works but then it does not recognize GBSW.

Normally I use: ./install.com gnu huge X86_64 IFORT little_e
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 06:56 PM
I am able to use:

update inbfrq -1 imgfrq -1 ihbfrq 0 -
atom switch cdie vdw vswitch cutnb 17 ctofnb 16 ctonnb 16

with GBSW and xxlarge. I guess this is enough for my Monte Carlo simulations since my protein and surface are rigid, so the list of neighbors does not need to be more than 1 A larger than the cutoff.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/27/11 06:59 PM
I'm still not sure it's going to help you. But it's always good to try! So, the version I built is pretty much useless to you because I got rid of all the implicit models.

What you need is in the following thread and links in it:

http://www.charmm.org/ubbthreads-7-5-5/ubbthreads.php?ubb=showflat&Number=25610

What I would try is adding -mcmodel=medium to the appropriate compiler options in /build/Makefile_gnu

it seems like CHARMM is trying to allocate a chunk of HEAP that's too large, because of the size of nonbonded list.
Posted By: lennart Re: Generalized Born and carbohydrates - 04/27/11 07:08 PM
Since simple truncation usually is silly this is not directly available in CHARMM; the specified settings are a trick to implement truncation.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 07:08 PM
Quote:
What I would try is adding -mcmodel=medium to the appropriate compiler options in /build/Makefile_gnu

-mcmodel=medium is already set per default when using gfortran. It still fails to link. This is from Makefile_gnu:

# gfortran (replacement for g77)
ifdef GFORTRAN
ifdef X86_64
FC = gfortran -DGNU -fdefault-integer-8 -frecord-marker=4 -mcmodel=medium
LD = gfortran -DGNU -fdefault-integer-8 -frecord-marker=4 -mcmodel=medium
else
FC = gfortran -DGNU
LD = gfortran -DGNU
endif
endif
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 07:24 PM
Originally Posted By: lennart
Since simple truncation usually is silly this is not directly available in CHARMM; the specified settings are a trick to implement truncation.

If the protein is rigid, how important is cutnb? Can it be set equal or just 1 A larger than the cutoff?

Thanks.
Posted By: rmv Re: Generalized Born and carbohydrates - 04/27/11 07:48 PM
The seg fault upon HEAP expansion can indicate a problem with the compilation; version 4.2.x or later of gfortran is highly recommended.

Reducing CUTNB can help with memory usage, at the expense of speed.
Posted By: lennart Re: Generalized Born and carbohydrates - 04/27/11 08:01 PM
Why would it make a difference that the protein is rigid? Other components still move, don't they? The most efficient setting of CUTNB (with INBFRQ -1 [IMGFRQ -1]) is easily determined for your particular conditions: Three runs of 1000 steps with CUTNB = CTOFNB +2, +4, +6 (and for parallel runs it is the elapsed time, not the CPU time that is important).

With INBFRQ -1 the choice of CUTNB is merely an efficiency issue. Some other program has (or by now hopefully had?) CUTNB=CTOFNB with INBRFQ >1, which means you HAVE to use a thermostat to hide the constant heating this creates.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 08:10 PM
Originally Posted By: rmv
The seg fault upon HEAP expansion can indicate a problem with the compilation; version 4.2.x or later of gfortran is highly recommended.

I am using 4.4 on one PC and 4.5 on another.

Originally Posted By: rmv

Reducing CUTNB can help with memory usage, at the expense of speed.

ok, I can live with that for the moment.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 08:14 PM
Originally Posted By: lennart
Why would it make a difference that the protein is rigid? Other components still move, don't they?

It's a Monte Carlo simulation. The system consists of a protein and a surface. Both are rigid. At each Monte Carlo move the protein is randomly rotated and translated. There is no dynamics involved.
Posted By: lennart Re: Generalized Born and carbohydrates - 04/27/11 08:47 PM
CUTNB obviously is irrelevant for MC (well, depending on the move set it might still be useful I guess).
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 09:34 PM
Also, I have noticed that using GBSW (with 17/16/16 cutoff) is only 2x slower than using the inaccurate RDIE method (with 14/12/10 cutoff). Also, GBSW is 25% faster than GBMV method II (same cutoff). Please, note that this might be completely different with dynamics simulations. I am running Monte Carlo for now and "all" I need is to evaluate the energy.

Thanks to all! This thread has been very helpful. Even if I am still not able to compile the huge version I have found a set of cutoff parameters which work well.
Posted By: lennart Re: Generalized Born and carbohydrates - 04/27/11 09:36 PM
If you need the huge version for implicit solvent (GB) calculations, it may well be that explicit solvent is more efficient....
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/27/11 09:45 PM
Originally Posted By: lennart
If you need the huge version for implicit solvent (GB) calculations, it may well be that explicit solvent is more efficient....

I'm not doing dynamics though. Monte Carlo is appropriate for the kind of sampling that I am looking for and I don't think I can perform Monte Carlo in explicit water.
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/28/11 01:13 AM
Originally Posted By: lennart
CUTNB obviously is irrelevant for MC (well, depending on the move set it might still be useful I guess).

By "irrelevant" do you mean that it could also be set equal to or even smaller than CTOFNB?

Thanks!
Posted By: lennart Re: Generalized Born and carbohydrates - 04/28/11 05:04 AM
No, that will not work.
Posted By: apredeus Re: Generalized Born and carbohydrates - 04/28/11 05:16 AM
I think CHARMM crashes if you have CUTNB < CTOFNB
Posted By: gianluca Re: Generalized Born and carbohydrates - 04/28/11 05:19 AM
I decided to test it. In my case, it does not crash but it produces nonsense smile
Posted By: Kenno Re: Generalized Born and carbohydrates - 04/29/11 02:39 AM
Originally Posted By: gianluca
I don't think I can perform Monte Carlo in explicit water.
Why not?
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/06/11 11:07 PM
Originally Posted By: Kenno
Originally Posted By: gianluca
I don't think I can perform Monte Carlo in explicit water.
Why not?

Don't you have to re-equilibrate the water molecules around the protein after each move? Or am I wrong? I would be happy for any advice. Thanks.
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/06/11 11:11 PM
UPDATE: I was able to compile the huge size of CHARMM by removing the options "-fdefault-integer-8 -frecord-marker" but still keeping "-mcmodel=medium".
Posted By: rmv Re: Generalized Born and carbohydrates - 05/06/11 11:29 PM
The "-frecord-marker=4" should have no effect on executable size, but may impact trajectory compatibility.
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/06/11 11:48 PM
Originally Posted By: rmv
The "-frecord-marker=4" should have no effect on executable size, but may impact trajectory compatibility.

I tried to leave "-frecord-marker=4" in, but it did not compile.
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/07/11 01:04 AM
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!
Posted By: Kenno Re: Generalized Born and carbohydrates - 05/07/11 01:33 AM
Originally Posted By: gianluca
Originally Posted By: Kenno
Originally Posted By: gianluca
I don't think I can perform Monte Carlo in explicit water.
Why not?

Don't you have to re-equilibrate the water molecules around the protein after each move? Or am I wrong? I would be happy for any advice. 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:
Originally Posted By: gianluca
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.
Originally Posted By: gianluca
The errors are bigger with larger Monte Carlo moves.
Your setup is unusual to me, but wouldn't all of that that be expected behavior?

Originally Posted By: gianluca
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.
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.

Originally Posted By: gianluca
I have played around with all possible parameters but can't figure it out.
Trying to gain insight into the problem by reading might in this case be more effective than a manual MC search of parameter space... grin
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/07/11 03:43 AM
Kenno,

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

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.
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/07/11 08:05 PM
Quote:
Trying to gain insight into the problem by reading might in this case be more effective than a manual MC search of parameter space...

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!
Posted By: rmv Re: Generalized Born and carbohydrates - 05/07/11 09:46 PM
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 ...
Posted By: Kenno Re: Generalized Born and carbohydrates - 05/10/11 12:45 AM
Originally Posted By: gianluca
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.
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":
Quote:
In the interest of computational efficiency, Monte Carlo calls specific
energy routines directly, rather than through the main ENERGY routine. As a
result, not all energy terms are supported. Those that are supported are
bonds, angles, Urey-Bradley, dihedrals, impropers, vdw, electrostatic,
image vdw, image electrostatic, QM/MM (MOPAC only), path integral, asp-EEF1,
asp-ACE/ACS, NOE constraints, and user (also note that the user must edit
both usersb.src and mcuser.src for the user energy to work correctly). All
non-bonded calculations can be either atom- or group-based. Terms not listed
above are not included in the present implementation.
Wait a second... I don't see GBMV on this list... looks like I made the right quip for the wrong reason. laugh

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. wink
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/10/11 02:23 AM
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
Posted By: Kenno Re: Generalized Born and carbohydrates - 05/10/11 03:50 AM
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?! grin
Posted By: gianluca Re: Generalized Born and carbohydrates - 05/10/11 04:16 AM
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. crazy )
© CHARMM forums