CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics Free energy between two graphenes with WHAM

Dear all,


I applied the restraints to two graphenes for calculation of the free energy as following:

**********************************************************************************************************

! Block to translate with x,y directions
mmfp
geo cylinder rcm -
xref 0.0 yref 0.0 zref 0.0 zdir 1.0 -
harm force 100 droff 1.0 sele segid GRAP .and. resi 1 end
end

mmfp
geo cylinder rcm -
xref 0.0 yref 0.0 zref 0.0 zdir -1.0 -
harm force 100 droff 1.0 sele segid GRAP .and. resi 2 end
end


! For reaction coordinates (Distance between two graphenes // Distance = 10~20 ang. with 21 windows)
mmfp
geo plane -
harm force 2 droff 1.0 -
xref 0.0 yerf 0.0 zref 10.00 -
xdir 0.0 ydir 0.0 zdir 1.0 -
sele segid GRAP .and. resi 1 .and. type C* end
end

mmfp
geo plane -
harm force 2 droff 1.0 -
xref 0.0 yerf 0.0 zref -10.00 -
xdir 0.0 ydir 0.0 zdir 1.0 -
sele segid GRAP .and. resi 2 .and. type C* end
end

********************************************************************************************************

I applied the plane restraints to each carbon atom because of fixing graphene's orientation.

After perform MD simulations, I obtained the free energy with WHAM software made by Prof. Grossfield (url: http://membrane.urmc.rochester.edu/?page_id=59 ).

When I ran the software, I set the force constant as 4.0 kcal/mol.

As a result, the free energy landscape shape looks well but scale of the energy difference is very low (0.5~1 kcal/mol, please find an attchment).

So, I wonder what is a problem.



Thank you to read my question.



P.S. two graphene have same structure, consist of 88 carbons.
And the size is 1.5x1.5 nm^2, performed NPT simulation with SATP condition.

It seems to be me that the cylindrical and planar MMFP restraints may conflict with each other, but it is difficult to tell from the script fragment posted.

Also, it should be noted that the DOMDEC code for faster MD does not support the use of MMFP restraints.

Thanks to reply.


First, I didn't use the DOMDEC codes.

Second, the results seem to be that the restraints are well working by checking the configurations over time.

In my script, I just called/set the force fields/the crystal structure/the system information/nonbonded part before the restraint part.
Also, after the restraint part, there is the dynamics part.
I made 21 independent systems (windows) with the different distance (10~20 ang.) between two graphenes with 0.5 ang. interval.

The cylindrical restraints are just for making graphene impossible to translate with x, y-directions,
which means that the reaction coordinate in this system is only the z-direction distance between two graphenes.

The reason why I use the plane restraints, is that normal vector of the graphenes must be perpendicular (or close to) to xy-plane for my study.


My problem is that The free energy landscape looks good, but the free energy difference between local minimum and other regions (just 0.5 kcal/mol).
Do I need to product the number of atoms 88 (88 atoms in 1 graphene) to the free energy landscape? (if right, then the free energy scale is too large.)


Thank you to read my question.