* GENERATED BY CHARMM-GUI (http://www.charmm-gui.org) v2.0 on Feb, 07. 2020. JOBID=8108355655
* INPUT FILES FOR PRODUCTION
* PRODUCTION COMMAND: $CHARMMEXEC CNT=1 < STEP7_PRODUCTION.INP > STEP7_PROUDCTION.OUT
*

DIMENS CHSIZE 3000000 MAXRES 3000000

! Read topology and parameter files
stream toppar.str

! Read the system information
stream step6.1_equilibration.str

! Read PSF
open read unit 10 card name step5_assembly.psf
read psf  unit 10 card xplor

! Read Coordinate
open read unit 10 card name step5_assembly.crd
read coor unit 10 card

!
! Previous box information
!

calc pcnt = @cnt - 1
bomlev -5
if pcnt .eq. 0 open read  unit 11 card name step6.6_equilibration.rst
if pcnt .gt. 0 open read  unit 11 card name step7_@pcnt.rst
read coor dynr curr unit 11
bomlev  0

calc A  = ?XTLA
calc B  = ?XTLB
calc C  = ?XTLC

!
! Image Setup
!
open read unit 10 card name crystal_image.str
CRYSTAL DEFINE @XTLtype @A @B @C @alpha @beta @gamma
CRYSTAL READ UNIT 10 CARD

!Image centering by residue
IMAGE BYRESID XCEN 0.0 YCEN 0.0 ZCEN @zcen sele resname TIP3 .or. segid MEMB end
IMAGE BYRESID XCEN 0.0 YCEN 0.0 ZCEN @zcen sele resname @posid .or. resname @negid end
define PROT sele ( segid PROA ) end
IMAGE BYSEGID XCEN 0.0 YCEN 0.0 ZCEN @zcen sele PROT end

!
! Nonbonded Options
!

system "python checkfft.py ?XTLA ?XTLB ?XTLC > checkfft.str"
stream checkfft.str

nbonds atom vatom vfswitch bycb -
       ctonnb 10.0 ctofnb 12.0 cutnb 16.0 cutim 16.0 -
       inbfrq -1 imgfrq -1 wmin 1.0 cdie eps 1.0 -
       ewald pmew fftx @fftx ffty @ffty fftz @fftz  kappa .34 spline order 6
energy

! To remove the draft of the protein
MMFP
GEO rcm cylinder -
    Xref 0.0 Yref 0.0 Zref 0.0 XDIR 0.0 YDIR 0.0 ZDIR 1.0 -
    harmonic FORCE 1.0 select .not. hydrogen .and. PROT end
END

!
! Production
!------------------------------------------------------------------------
!
! ** Note: change "nstep" to reduce the number of dynamics steps
!

set nstep = 500000

! estimate Pmass from SYSmass (total system mass)
! [there could be problems with exreme values, such as  Pmass << SYSmass or Pmass >> SYSmass
scalar mass stat
calc Pmass = int ( ?stot  /  50.0 )

set temp = 323.15
shake bonh param fast
calc pcnt = @cnt - 1
if pcnt .eq. 0 open read  unit 11 card name step6.6_equilibration.rst
if pcnt .gt. 0 open read  unit 11 card name step7_@pcnt.rst
open write unit 12 card name step7_@cnt.rst
open write unit 13 file name step7_@cnt.dcd

DYNA CPT leap restart nstep @nstep timestep 0.002 echeck -1 -
     nprint 1000 iprfrq 5000 ntrfrq 5000 -
     iunrea 11 iunwri 12 iuncrd 13 iunvel -1 kunit -1 -
     nsavc 50000 nsavv 0 -
     PCONS pint pref 1.0 pmass @Pmass pgamma 20.0 - 
     HOOVER reft @temp tmass 2000.0 tbath @temp firstt @temp

!
! write some information
!
open write unit 10 card name step7_@cnt.pdb
write coor unit 10 pdb
close unit 10

open write unit 10 card name step7_@cnt.crd
write coor unit 10 card
close unit 10
stop
