remarks  h2orefine.inp
remarks  Michael Nilges, FEB-1994, EMBL
remarks  Jens Linge, AUG-1998, EMBL


evaluate ($filename=$file)
structure reset @@$structure end

coor @@$filename
vector do (refx = x) (all)
vector do (refy = y) (all)
vector do (refz = z) (all)

vector do (segid = "PROT") (segid "    ")

set message on echo on end

@generate_water.xplor

vector do (segid = "    ") (segid "PROT")


@@read_data_test.xplor

flags exclude * end
@@flags_h2o.xplor
flags include elec end


noe scale * 50 end
evaluate ($kcdih = 200)
restraints dihedral 
      scale=$kcdih                      
end 


@setup_swap.xplor
parameter
   improper (store1) (store1) (all) (all) 0.0 TOKEN TOKEN
   improper (all) (all) (store1) (store1) 0.0 TOKEN TOKEN
end


! since we do not use SHAKe, we increase the water bond angle energy constant
parameter 
   angle (resn tip3) (resn tip3) (resn tip3) 500 TOKEN 
end

! reduce angle force constant for some atoms
evaluate ($kangle = 50)
evaluate ($kimpro = 5)
vector ident (store1) (store1 or resn tip3)
parameter
   angle    (not resn tip3)(not resn tip3)(not resn tip3) $kangle  TOKEN
   improper (not store1)(not store1)(not store1)(not store1) $kimpro  TOKEN TOKEN
end

! in case we use planar restraints, we can switch off all impropers.
!parameter
!   improper (all) (all) (all) (all) 0.0 TOKEN TOKEN
!end
!evaluate ($kplanar = 5)
!@planarf.xplor

constraints fix (not segid w*) end
minimize powell nstep=20 drop=100 end

constraints fix (segid none) end
restraints harmonic 
  reference = main
  exponent = 2
end
vector do (harm = 0)  (all)
vector do (harm = 10) (name ca)

minimize powell nstep=20 drop=10 end
minimize powell nstep=20 drop=10 end

vector do (mass = 100) (all)
vector do (fbeta = 0) (all)
vector do (fbeta = 20. {1/ps} ) ( not name h* )                

for $bath in (100 200 300 400 500) loop heat
   dynamics langevin
     nstep=50 timest=0.005 {ps}      
     iasvel=maxwell  firsttemperature = $bath 
     tbath=$bath  tcoupling = true 
     nprint=50 iprfrq=50                ! statistics output
   end 
end loop heat

vector do (harm = 0)  (all)

dynamics langevin
     nstep=2500 timest=0.005 {ps}      
     iasvel=maxwell  firsttemperature = $bath 
     tbath=$bath  tcoupling = true 
     nprint=50 iprfrq=50                ! statistics output
end 

for $bath in (500 400 300 200 100) loop cool
   evaluate ($kangle = min(500,$kangle * 2))
   evaluate ($kimpro = min(500,$kimpro * 3))
   parameter
      angle    (not resn tip3)(not resn tip3)(not resn tip3) $kangle  TOKEN
      improper (not store1)(not store1)(not store1)(not store1) $kimpro  TOKEN TOKEN
   end
!  evaluate ($kplanar = 10) @SEQUENCE:planarf.xplor
   dynamics langevin
     nstep=200 timest=0.005 {ps}      
     iasvel=maxwell  firsttemperature = $bath
     tbath=$bath  tcoupling = true                        
     nprint=50 iprfrq=50                ! statistics output
   end 
end loop cool


mini powell nstep 100 end