#============================================================*
#               test lennard-jones potential                 *
#                                                            *
# use sidechain atoms.                                       *
#                                                            *
# specify a single list of potentials.                       *
#                                                            *
# test iq-iq repulsice potential.                            *
#============================================================*

system echo = off
system verbose = false


# set units
#----------
units
units scale  name = mass   uname = reducedMass   value = 10000.0

# open pdb file
#--------------
database open name =  neck_xform  format = pdb file =  neck_xform.pdb

# read molecule
#--------------
molecule read  database =  neck_xform  name = m1  model = 1  type = protein

# create domains
domain create hc = C
domain hc  backbone  color = [ 0.6 0.6 0.6 ] width = 2 show
domain create iq1 = AC[760-788]
domain iq1 backbone color = [ 1 0 0 ] width = 2 show
domain create iq2 = BC[788-814]
domain iq2 backbone color = [ 0 1 0 ] width = 2 show
domain create iq3 = DC[814-836]
domain iq3 backbone color = [ 1 0.6 0 ] width = 2 show

#--------------------------------------------------------------#
#                      create bodies                           #
#--------------------------------------------------------------#
# 1st group
body create  name = ground    type = ground
body create  name = iq1_body  type = rigid   domain = iq1
body create  name = iq2_body  type = rigid   domain = iq2
body create  name = iq3_body  type = rigid   domain = iq3


#--------------------------------------------------------------#
#                      create joints                           #
#--------------------------------------------------------------#
# 1st group
joint create  name = jnt1   type = weld    domain = hc    res = C[760]     use = CA
joint jnt1 bodies = { ground iq1_body }
joint jnt1 color = [ 1 0 1 ] msize = 0.1 shading = color show = true 

joint create  name = jnt2   type = ball    domain = hc    res = C[788]     use = CA
joint jnt2 bodies = { iq1_body iq2_body }
joint jnt2 color = [ 1 0 1 ] msize = 0.1 shading = color show = true

joint create  name = jnt3   type = ball    domain = hc    res = C[814]     use = CA
joint jnt3 bodies = { iq2_body iq3_body }
joint jnt3 color = [ 1 0 1 ] msize = 0.1 shading = color show = true 

#--------------------------------------------------------------#
#                 create a rigid simulation                    #
#--------------------------------------------------------------#

simulation create name = rsim  type = rigid

# add bodies
#-----------
simulation rsim add body = ground
simulation rsim add body = iq1_body
simulation rsim add body = iq2_body
simulation rsim add body = iq3_body

body iq1_body  damping = 20.0
body iq2_body  damping = 20.0
body iq3_body  damping = 20.0


# add joints
#-----------
simulation rsim add joint = jnt1
simulation rsim add joint = jnt2
simulation rsim add joint = jnt3


# set solver
#-----------
simulation rsim solver = ode
simulation rsim solver finite_rotation = on


# add forces
#-----------

force create name = force1           \
             type = explicit         \
             point = [ 0 0 0 ]       \
             direction = [ 8 -8 0 ]  \
             scale = 1.0

simulation rsim add force = force1    \
                    body = iq3_body   \
                    time = { 0 100 }


# add interaction potentials
#---------------------------

variable  show = true  

#domain iq1 define region = iq1_rgn  residues = C[760-787]  no_mainchain = true 
domain iq1 define region = iq1_rgn  residues = C[760-787]A[1-148]  no_mainchain = true 
#domain iq1 define region = iq1_rgn  residues = A[1-148]  no_mainchain = true 

#domain iq2 define region = iq2_rgn  residues = C[788-813] no_mainchain = true  
domain iq2 define region = iq2_rgn  residues = C[788-813]B[1-148]  no_mainchain = true  
#domain iq2 define region = iq2_rgn  residues = B[1-148]  no_mainchain = true  

domain iq3 define region = iq3_rgn  residues = C[814-836]D[1-148]  no_mainchain = true  
#domain iq3 define region = iq3_rgn  residues = D[1-148]  no_mainchain = true  

body iq1_body add potential = iq1_mm                                         \
     type = molecular-mechanics  region = iq1_rgn  geometry = points        \
     display_spheres = true  display = solid color = [ 1 0 0 ]  show = ${show}

body iq2_body add potential = iq2_mm                                         \
     type = molecular-mechanics  region = iq2_rgn  geometry = points        \
     display_spheres = true  display = solid color = [ 0 1 0 ]  show = ${show}

body iq3_body add potential = iq3_mm                                        \
     type = molecular-mechanics  region = iq3_rgn  geometry = points       \
     display_spheres = true  display = solid  color = [ 1 0.7 0 ]   show = ${show}            
#----------#

simulation rsim add interaction = mm1               \
                    list = { iq1_mm iq2_mm iq3_mm } \
                    time = { 0.0 10.0 }


# initialize simulation
#----------------------
simulation rsim initialize
simulation rsim momentum = on
simulation rsim damping = on

simulation rsim time step = 0.002
simulation  state  save_frequency = 20
#simulation  rsim  write  energy = true  contact = true  file_name = rsim

graphics  center  point = [ 1.74683 3.46701 2.54842 ]