remarks   Asn MD with GB/HCT
remarks   Tom Simonson, May, 2002. 


!--------------------------------------------------------------
! the initial time is the ONLY parameter that has to be changed
! $inittime is set to the last time of the previous run

evaluate ($inittime = 0 {ps})
evaluate ($finaltime = $inittime + 100)  {100 ps equilibration}
set seed = 123456 end


topology  
  @TOPPAR:topamber.inp 
  @TOPPAR:amberpatches.pro
end

parameter @TOPPAR:paramber.gb.inp end


segment
name="ASN1"
molecule name=ASN number=1 end
end
patch NASN refe=nil=(resid 1) end
patch CASN refe=nil=(resid 1) end


vector do (charge = charge + .2306) (name ht%)
vector show sum (charge) (all)



segment
name="ASN1"
molecule name=ASN number=1 end
end
patch NTER refe="+"=(resid 1) end
patch CTER refe="-"=(resid 1) end

vector show sum (charge) (all)



parameter 
nbonds 
    atom cdie trunc 
    e14fac=0.8333333                     ! use this to reproduce amber elec
    cutnb 500. ctonnb 480. ctofnb 490.   ! essentially no cutoff
    tolerance=100.                       ! only build the nonbonded list once
    nbxmod 5 vswitch
    wmin=1.0
end 
end

parameters nbonds
    EPS=1. WEPS=80. offset=0.09  GBHCT  ! GB parameters
end 
end




! gauche minus conformer
coor ! @data/asn.-80.10.pdb
REMARK  CHI1=-80, CHI2=10; ENERGY = -3.92838                                          
REMARK   DATE:     3/ 6/98     14:38:33      CREATED BY USER: simonson                
ATOM      1  N   ASN     1      26.522 -15.164  -2.651  0.00  0.00      ASN1
ATOM      2  HT1 ASN     1      25.821 -14.402  -2.738  0.00  0.00      ASN1
ATOM      3  HT2 ASN     1      26.992 -15.315  -3.566  0.00  0.00      ASN1
ATOM      4  HT3 ASN     1      27.228 -14.908  -1.932  0.00  0.00      ASN1
ATOM      5  CA  ASN     1      25.834 -16.419  -2.236  0.00  0.00      ASN1
ATOM      6  HA  ASN     1      25.080 -16.641  -2.979  0.00  0.00      ASN1
ATOM      7  CB  ASN     1      26.844 -17.593  -2.073  0.00  0.00      ASN1
ATOM      8  HB1 ASN     1      26.370 -18.398  -1.469  0.00  0.00      ASN1
ATOM      9  HB2 ASN     1      27.757 -17.253  -1.543  0.00  0.00      ASN1
ATOM     10  CG  ASN     1      27.223 -18.202  -3.434  0.00  0.00      ASN1
ATOM     11  OD1 ASN     1      26.897 -17.663  -4.497  0.00  0.00      ASN1
ATOM     12  ND2 ASN     1      27.947 -19.361  -3.377  0.00  0.00      ASN1
ATOM     13 HD21 ASN     1      28.231 -19.812  -4.222  0.00  0.00      ASN1
ATOM     14 HD22 ASN     1      28.188 -19.756  -2.491  0.00  0.00      ASN1
ATOM     15  C   ASN     1      25.143 -16.177  -0.888  0.00  0.00      ASN1
ATOM     16  OT1 ASN     1      23.884 -16.213  -0.858  0.00  0.00      ASN1
ATOM     17  OT2 ASN     1      25.860 -15.961   0.127  0.00  0.00      ASN1
END

write coor output=asn.-80.10.pdb end


@TOPPAR:volumes.amber
vector do (rmsd = rmsd * 0.9) (all)  ! reduce volumes by 10%
flags include gbse gbin end
parameter reduce selection=(all) overwrite=true mode=average end end

energy end


!-----------------------------
! shake TIP3p's and ALL bonds 
!shake
!    ncon=6000
!    reference=parameter bonds ( not resname TIP3 )  ( all )
!    tolerance=1.0e-08
!end


! weakly restrain Calpha position
vector do (harmonic = 0.5) (name ca)
coor disp=reference @asn.-80.10.pdb
constraints harmonic end
flags include harmonic end



!-----------------
! Now run dynamics 


!------------------------------------------
! $time is set to last time of previous run

evaluate ($time = $inittime)


!---------------------------------------------------------------
! first, read "restart" files or initial coords and assign velos
  


if ($time > 0) then

  coor disp=main @@$fincrd
  coor disp=comp @@$finvel

  vector do (vx = xcomp) (attribute xcomp < 999.0) 
  vector do (vy = ycomp) (attribute xcomp < 999.0)
  vector do (vz = zcomp) (attribute xcomp < 999.0)
  vector do (vx = maxwell(250)) (attribute xcomp > 999.)
  vector do (vy = maxwell(250)) (attribute xcomp > 999.)
  vector do (vz = maxwell(250)) (attribute xcomp > 999.)

else
  
  
  vector do (vx = maxwell(250)) (all)
  vector do (vy = maxwell(250)) (all)
  vector do (vz = maxwell(250)) (all)

end if



!-------------------------------------
! relax initial bad contacts
if ($time = 0.) then
minimize powell nstep=10 drop=10 end
minimize powell nstep=10 drop=10 end
end if



!--------------------------
! generate output filenames
  
!evaluate ($traj = "output/dyn" + encode($time) + ".crd")
!evaluate ($velo= "output/dyn" + encode($time) + ".vel")
!evaluate ($restart ="output/dyn"+encode($time) + ".rst")


!--------------
! let it go

display  CPU= $CPU sec.


dynamics verlet
     nstep=50 timest=0.001 {ps}     ! 100 ps dynamics
     iasvel=current                     ! current velocities 
     nprint=25 iprfrq=250              ! statistics output
     ascii=true
     !traj=$traj                         ! coord trajectory
     nsavc=250
     vascii=true voffset=800
     !velo=$velo                         ! velocity trajectory
     nsavv=10000
     !save=$restart
     isvfrq=100000                       ! restart file
end 




stop