# This example looks complicated because I split the
# hexadecane molecule into individual CH2 and CH3 monomers.
#
# I defined it this way so that you can easily modify 
# it to change the length of the alkane chain.


import "gaff.lt"      # load the "GAFF" force-field information
import "ch2group.lt"  # load the definition of the "CH2" object
import "ch3group.lt"  # load the definition of the "CH3" object



Hexadecane inherits GAFF {


  # Create an array of 16 "CH2" objects distributed along the X axis

  monomers = new CH2.move(0,0.4431163,0) [16].rot(180,1,0,0).move(1.2533223,0,0)

  #   Each CH2 monomer is initial moved in the +Y direction by 0.43116
  #   angstroms. Then it is rotated 180 degrees with respect to the
  #   previous monomer, and moved 1.2533223 Angstroms down the X axis.

  # ---- Now, modify the ends: ---
  # Delete the CH2 groups at the beginning and end, and replace them with CH3.
  # (Note: Alternately, instead of deleting the CH2 groups at each end, you
  #  could modify them by adding an extra hydrogen atom to those carbons.)

  delete monomers[0]
  delete monomers[15]
  endcap1 = new CH3
  endcap2 = new CH3

  # Move the CH3 groups to the correct location at either end of the chain:

  endcap1.move(0,0.4431163,0)
  endcap2.move(0,0.4431163,0).rot(180,0,0,1).move(18.7998345,0,0)

  # Note: 18.7998345 = (16-1) * 1.2533223


  # Now add a list of bonds connecting the carbon atoms together:

  write('Data Bond List') {
    $bond:b1  $atom:endcap1/C $atom:monomers[1]/C
    $bond:b2  $atom:monomers[1]/C  $atom:monomers[2]/C
    $bond:b3  $atom:monomers[2]/C  $atom:monomers[3]/C
    $bond:b4  $atom:monomers[3]/C  $atom:monomers[4]/C
    $bond:b5  $atom:monomers[4]/C  $atom:monomers[5]/C
    $bond:b6  $atom:monomers[5]/C  $atom:monomers[6]/C
    $bond:b7  $atom:monomers[6]/C  $atom:monomers[7]/C
    $bond:b8  $atom:monomers[7]/C   $atom:monomers[8]/C
    $bond:b9   $atom:monomers[8]/C  $atom:monomers[9]/C
    $bond:b10  $atom:monomers[9]/C  $atom:monomers[10]/C
    $bond:b11  $atom:monomers[10]/C  $atom:monomers[11]/C
    $bond:b12  $atom:monomers[11]/C  $atom:monomers[12]/C
    $bond:b13  $atom:monomers[12]/C  $atom:monomers[13]/C
    $bond:b14  $atom:monomers[13]/C  $atom:monomers[14]/C
    $bond:b15  $atom:monomers[14]/C  $atom:endcap2/C
  }

                       # OPTIONAL:
  create_var { $mol }  # Create a molecule ID number.  This number will
                       # be shared by all of the atoms in this polymer.
                       # In ch2group.lt, "$mol:..." refers to this number.

} # Hexadecane










######### (scratchwork calculations for the atomic coordinates) #########
#
#   1.2533223 = DeltaXc   = how far each CH2 group is shifted along 
#                           the X axis (in Angstoms).
#   0.4431163 = DeltaYc/2 = lateral displacement of carbons away
#                           from the central axis.  (See below.)
#
# Lcc = 1.5350 # length of the C-C bond (Sp3)
# Lch = 1.0930 # length of the C-H bond
# theta=2*atan(sqrt(2)) # ~= 109.5 degrees = tetrahedronal angle (C-C-C angle)
# DeltaXc = Lcc*sin(theta/2)  # = 1.2533222517240594
# DeltaYc = Lcc*cos(theta/2)  # = 0.8862326632060754
#                             # 0.5*DeltaYc = 0.4431163316030377
# DeltaZh = Lch*sin(theta/2)  # = 0.8924307629540046
# DeltaYh = Lch*cos(theta/2)  # = 0.6310438442242609