remarks dsele.inp
remarks test double selection facilities
remarks Thomas Simonson
remarks ================

topology  @CNS_TOPPAR:protein.top end
parameter @CNS_TOPPAR:protein.param end
parameter
  nbonds
     atom cdie shift eps=1.0
     cutnb=100.0
     nbxmod=5 vswitch 
  end
end

segment
  name=tryp
  chain
     @CNS_TOPPAR:protein.link
     sequence
        TRP
     end
  end
end

coor
ATOM      1  HT1 TRP     1       0.101   0.574  -1.385  1.00  0.00      TRYP
ATOM      2  HT2 TRP     1       1.711   0.436  -1.904  1.00  0.00      TRYP
ATOM      3  N   TRP     1       1.067   0.302  -1.101  1.00  0.00      TRYP
ATOM      4  HT3 TRP     1       1.080  -0.689  -0.791  1.00  0.00      TRYP
ATOM      5  CA  TRP     1       1.472   1.192   0.027  1.00  0.00      TRYP
ATOM      6  CB  TRP     1       2.696   2.026  -0.348  1.00  0.00      TRYP
ATOM      7  CG  TRP     1       3.922   1.169  -0.679  1.00  0.00      TRYP
ATOM      8  CD2 TRP     1       4.295   0.679  -1.890  1.00  0.00      TRYP
ATOM      9  CE2 TRP     1       5.424  -0.010  -1.601  1.00  0.00      TRYP
ATOM     10  CE3 TRP     1       3.769   0.774  -3.155  1.00  0.00      TRYP
ATOM     11  CD1 TRP     1       4.813   0.775   0.290  1.00  0.00      TRYP
ATOM     12  NE1 TRP     1       5.727   0.046  -0.320  1.00  0.00      TRYP
ATOM     13  HE1 TRP     1       6.094  -0.743   0.127  1.00  0.00      TRYP
ATOM     14  CZ2 TRP     1       6.138  -0.677  -2.563  1.00  0.00      TRYP
ATOM     15  CZ3 TRP     1       4.472   0.108  -4.148  1.00  0.00      TRYP
ATOM     16  CH2 TRP     1       5.636  -0.605  -3.856  1.00  0.00      TRYP
ATOM     17  C   TRP     1       0.280   2.113   0.309  1.00  0.00      TRYP
ATOM     18  OT1 TRP     1      -0.402   2.481  -0.636  1.00  0.00      TRYP
ATOM     19  OT2 TRP     1      -0.106   2.192   1.459  1.00  0.00      TRYP
END

energy end
pick bond (name CA) (name CB) energy
pick angle (name CA) (name CB) (name CG) energy
pick angle (name CA) (name CB) (name CG) geometry
print threshold=3. angles

{* subdivide molecule into two parts *}
{=====================================}

igroup
 interaction       (name HT* or name N or name CA or name C or name OT*)
                   (prev)
 weights bond 0.10 angle 0.20 dihe 0.30 impr 0.40 vdw 0.50 elec 0.60 end
 interaction       (name HT* or name N or name CA or name C or name OT*)
                   (prev)
 weights  bond 0.90 angle 0.80 dihe 0.70 impr 0.60 vdw 0.50 elec 0.40 end
 interaction       (name HT* or name N or name CA or name C or name OT*)
                   (not (prev))
 interaction       (not (name HT* or name N or name CA or name C or name OT*))
                   (prev)
 ?
end

energy end
pick bond (name CA) (name CB) energy
pick angle (name CA) (name CB) (name CG) energy
pick angle (name CA) (name CB) (name CG) geometry
print threshold=3. angles
show element (DX) (all)
show element (DY) (all)
show element (DZ) (all)


{* Create a perturbation by placing an additional charge on HT1 *}
{================================================================}

igroup 
   interaction (name HT1) (not name HT1)
   interaction (not name HT1) (not name HT1) 
end

set seed=271828. end
do (fbeta=100.) ( all ) 
do (vx=maxwell(298./2000)) ( all )
do (vy=maxwell(298./2000)) ( all )
do (vz=maxwell(298./2000)) ( all )

dyna tors
   nstep=10
   cmremove=false
   topology
      maxchn=300
      kdihmax=10.
   end
   timestep=0.003
   tcoupling=true   temperature=298.
   trajectory=coor.dat nsavc=1
end


{now calculate the free energy change "by hand" through the command language}
igroup interaction (name HT1) (not name HT1) end
energy accum=reset end
flags 
  exclude *
  include elec
end

{initializations for post-averaging of exp(-V/kT)}
read trajectory begin=1 skip=1 stop=10
     input=coor.dat
end
energy end
evaluate ( $N=1)
evaluate ( $V =-$ELEC/(298.*$KBOLTZ) )
evaluate ( $E = exp($V) )

{loop for post-averaging}

while ($N < 10) loop main
   read trajectory next end
   energy end
   evaluate ( $V =-$ELEC/(298.*$KBOLTZ) )
   evaluate ( $E = $E + exp($V) )
   evaluate ( $N = $N+1.)
end loop main
evaluate ( $E = $E/$N)
evaluate ( $A = LOG($E))
evaluate ( $A =-298.*$KBOLTZ*$A)
display $A  { free energy change }

igroup ? end

stop