Active Threads | Active Posts | Unanswered Today | Since Yesterday | This Week
Energy terms, Constraints, Restraints, and Solvation Jump to new posts
Re: syntax for simple restraint to box center lennart 10/23/20 05:15 PM
coor hbond can handle simple periodic boundary conditions (corman.info).
How good are the FF parameters for hydronium (a very rare species)?
4 70 Read More
Minimization, Normal modes, Monte Carlo,... Jump to new posts
GCMC simulation of crystal water Allen_123 10/23/20 10:40 AM
Hello Charmm users,

I'm currently using Grand canonical monte carlo ensemble to simulate crystal water by inserting and deleting water atoms inside a crystal. The crystal.crd contains position of crystal atoms as well as H2O. crystal atoms are fixed and set inactive. I got an error of
MC E ENR: Eval# ENERgy Delta-E GRMS
MC E INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
MC E EXTERN: VDWaals ELEC HBONds ASP USER
MC E IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
MC E EWALD: EWKSum EWSElf EWEXcl EWQCor EWUTil
MC E PRESS: VIRE VIRI PRESSE PRESSI VOLUme
MC E LRCor: EVDW VIRI
---------- --------- --------- --------- --------- ---------
MC E> 0 28693807.808 0.000 39991243.434
MC E INTERN> 39.71400 71.43027 9.33503 125.67478 2.99950
MC E EXTERN> 157.74161 -227.80326 0.00000 0.00000 0.00000
MC E IMAGES> 28694837.626 -689.860 0.000 0.000 0.000
MC E EWALD> 43.99046 -3024.23481 2466.52214 0.00000 0.00000
MC E PRESS> -0.52165E+10 0.11479E+09 0.13065E+12 0.28750E+10 0.27377E+04
MC E LRCor> -5.32850 -10.65528
---------- --------- --------- --------- --------- ---------

***** LEVEL 2 WARNING FROM *****
***** REACHED GCMC LIMIT
******************************************
BOMLEV ( 0) IS NOT REACHED. WRNLEV IS 5

My scripts are like following:

*gcmc test simulation
*
!--- Read parameters, psf and coordinates
open read card name "top_all36_prot.rtf" unit 25
read rtf card unit 25
close unit 25
open read card name "par_all36m_prot.prm" unit 25
read parameter card unit 25 flex
close unit 25
stream toppar_water_ions.str

set cutnb 16.0 ! cutnb
set ctonnb 10.0 ! ctonnb
set ctofnb 12.0 ! ctofnb
set eatom atom
set etrun switch
set vatom vatom
set vtrun vswitch


OPEN UNIT 1 FORM READ NAME unit.psf
READ PSF CARD UNIT 1
CLOSe UNIT 1

! Read coordinates of starting structure
OPEN READ CARD UNIT 1 NAME crystal.crd
READ COOR CARD UNIT 1
CLOSE UNIT 1

!Set the active and blocking atoms
SCALar GCMC SET 0.0 SELEct IRES 1:12 END ! ires 1 to 12 is crystal
SCALar GCMC SET 1.0 SELEct IRES 13:42 END ! ires 13 to 42 is water
SCALar GCBLocker SET 1.0 SELEct all end
crystal define HEXAgonal 24.0709 24.0709 5.4560 90.00 90.00 120.00
crystal build noper 0 cutoff 30
update inbfrq -1 imgfrq -1 ihbfrq 0 imall -
ewald pmewald kappa 0.34 fftx 32 ffty 32 fftz 32 order 6 lrc -
@eatom @etrun @vatom @vtrun cutnb @cutnb ctonnb @ctonnb ctofnb @ctofnb -
cutim 999
ENERGY

! Create the MC move set
! Rigid translations
MOVE ADD MVTP RTRN BYREsidue WEIGht 1.0 DMAX 0.25 -
SELE (TYPE OH2) END LABEl DISP
! Rigid rotations
MOVE ADD MVTP RROT BYREsidue WEIGht 1.0 DMAX 30.0 -
SELE (TYPE OH2) END LABEl ROTA
! Insertion and deletion
MOVE ADD MVTP GCMC WEIGht 1.0 SELE (TYPE OH2) END LABEl GCMC

! Link the GCMC moves to the rotations and displacements to avoid moving
! inactive molecules.
MOVE LINK GCMC LAB1 GCMC LAB2 DISP
MOVE LINK GCMC LAB1 GCMC LAB2 ROTA

! Link the translations and rotations to each other for greater efficiency
MOVE LINK LAB1 DISP LAB2 ROTA

OPEN WRITE UNFOrmatted UNIT 32 NAME gcmc.dcd
! Do 1000000 steps at 298 K, writing energy 1000 steps.
! Grid-based insertions with 10 attempted orientations are used.
MC NSTEp 1000000 TEMPerature 298.0 ISEEd 302430 -
INBFrq -1 IECHeck 1000 IMGFrq -1 IUNC 32 NSAVc 1000 imall -
MUEX -5.8 DENS 0.03342 -
RGRId 0.25 GCCUt 2.5 NOTB 10 -
XMIN -12 YMIN -12 ZMIN -2.7 -
XMAX 12 YMAX 12 ZMAX 2.7

DELEte ATOM SELEct .BYRES. PROP X .GT. 9998.0 END
open write card unit 44 name gcmc.crd
write coor card unit 44 sele gcmc end


Could any one tell where the issue is?

Thank you.
0 8 Read More
Molecular Dynamics Jump to new posts
Re: Free energy between two graphenes with WHAM trying_MD 10/23/20 06:38 AM
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.
2 24 Read More
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.006s Queries: 6 (0.003s) Memory: 0.8127 MB (Peak: 0.8420 MB) Data Comp: Off Server Time: 2020-10-24 03:55:28 UTC
Valid HTML 5 and Valid CSS