PForceField {


  # LAMMPS supports a large number of "styles" (ie. equations for calculating 
  # forces between particles).  At some point, we must eventually select the
  # formulas we want to use.  This can be done anywhere, but we might as
  # well specify that now.  Later on we will specify the parameters
  # which go into these equations.

  write_once("In Init") {
    # -- Styles used in "ForceField" --
    # -- (Changing these styles will change the formulas above) --
    units           real
    atom_style      full
    bond_style      class2
    angle_style     quartic
    dihedral_style  opls
    pair_style      buck 
  }


  # There are 4 atom types: "C" "H" "Si" "O"
  write_once("Data Masses") {
   @atom:Si   28.085
   @atom:C    12.01
   @atom:H     1.007
   @atom:O    15.999
  }
  # ---- 2-body (non-bonded) interactions: ----
  #   U(r) = Aexp[-r/B]-C/r^6
  #   (for details see http://lammps.sandia.gov/doc/pair_lj.html)
  #        atom-type   atom-type   epsilon sigma

  write_once("In Settings") {
    # Pairwise (non-bonded) interactions:
    #           atomType1 atomType2      A          B         C  
    pair_coeff   @atom:Si @atom:Si      50326      0.3513    3085.3
    pair_coeff   @atom:Si @atom:O       61781      0.2895    1109.4
    pair_coeff   @atom:Si @atom:C       24753      0.3369    1406.1
    pair_coeff   @atom:O  @atom:O       75844      0.2461     398.9
	pair_coeff   @atom:O  @atom:C       33702      0.2795     505.6
    pair_coeff   @atom:O  @atom:H       14177      0.2562     104.5
	pair_coeff   @atom:C  @atom:C       14967      0.3236     640.8
    pair_coeff   @atom:C  @atom:H        6300      0.2928     132.5
	pair_coeff   @atom:H  @atom:H        2650      0.2673      27.4

	
    # (Interactions between different atoms are determined by mixing rules.)
  }

  # ---- 2-body (bonded) interactions: ----
  #
  #   Ubond(r) = K2(r-r0)^2+K3(r-r0)^3+K4(r-r0)^4
  #   (for details see http://lammps.sandia.gov/doc/bond_harmonic.html)
  #
  write_once("In Settings") {
    #             bond-type        atomType1   atomType2        r0       K2       K3       K4 
    bond_coeff  @bond:Sidechain1    @atom:C     @atom:H         1.092     328      0        0        
    bond_coeff  @bond:Sidechain2    @atom:Si    @atom:C         1.878     190     -279      308
	bond_coeff  @bond:Backbone      @atom:Si    @atom:O         1.651     350     -517      674
  }

  # ---- 3-body angle (hinge) interactions: ----

  # Rules for determining 3-body (angle) interactions by atom & bond type:
  # angle-type     atomType1 atomType2 atomType3  bondType1 bondType2

  write_once("Data Angles By Type") {
    @angle:Sidechain4  @atom:H   @atom:C   @atom:H    @bond:*   @bond:*
	@angle:Sidechain3  @atom:Si  @atom:C   @atom:H    @bond:*   @bond:*
	@angle:Sidechain2  @atom:C   @atom:Si  @atom:C    @bond:*   @bond:*
	@angle:Sidechain1  @atom:O   @atom:Si  @atom:C    @bond:*   @bond:*
	@angle:Backbone2   @atom:O   @atom:Si  @atom:O    @bond:*   @bond:*
	@angle:Backbone1   @atom:Si  @atom:O   @atom:Si   @bond:*   @bond:*
    
  }

  # Force-field parameters for 3-body (angle) interactions:
  #
  #   Uangle(theta) = k2*(theta-theta0)^2+k3*(theta-theta0)^3+k4*(theta-theta0)^4
  #   (for details see https://docs.lammps.org/angle_quartic.html)
  #
  write_once("In Settings") {
    #             angle-type         theta0       K2         K3        K4
    angle_coeff @angle:Sidechain4    107.77      38.500      0         0 
    angle_coeff @angle:Sidechain3    111.09     -13.950      0         0
	angle_coeff @angle:Sidechain2    112.44      36.210     -20.39     20.020
    angle_coeff @angle:Sidechain1    109.82      23.022     -31.399    24.981
	angle_coeff @angle:Backbone2     105.56      91.835      0         0
    angle_coeff @angle:Backbone1     137.63      10.305     -18.101    10.100
  }

  # ---- 4-body dihedral interactions ----

  # 4-body interactions in this example are listed by atomType
  # Rules for determining 4-body (dihedral) interactions by atom & bond type:
  write_once("Data Dihedrals By Type") {
    # dihedralType atmType1 atmType2 atmType3 atmType4 bondType1 bType2 bType3
    @dihedral:Backbone @atom:Si @atom:O @atom:Si @atom:O  @bond:* @bond:* @bond:*
    @dihedral:Sidechain3 @atom:Si @atom:O @atom:Si @atom:C  @bond:* @bond:* @bond:*
	@dihedral:Sidechain2 @atom:O  @atom:Si @atom:C @atom:H  @bond:* @bond:* @bond:*
	@dihedral:Sidechain1 @atom:C  @atom:Si @atom:C @atom:H   @bond:* @bond:* @bond:*
  }

  # The forumula used is:
  #
  # Udihedral(phi) = (K1/2)*(1+cos(phi))  +  (K2/2)*(1+cos(2*phi)) +
  #             ...  (K3/2)*(1+cos(3*phi)) + (K4/2)*(1+cos(4*phi))
  #     The d parameter is in degrees, K is in kcal/mol/rad^2.
  #     (for details, see http://lammps.sandia.gov/doc/dihedral_charmm.html)
  #
  # The corresponding command is 
  # dihedral_coeff dihedralType               K1       K2        K3      K4

  write_once("In Settings") {
    dihedral_coeff @dihedral:Backbone       -0.41    0.27319    0.21968   0
    dihedral_coeff @dihedral:Sidechain3      0       0         -0.05706   0
	dihedral_coeff @dihedral:Sidechain2      0       0         -0.15      0
    dihedral_coeff @dihedral:Sidechain1      0       0         -0.15      0
  }

} # "ForceField"