#######################################################################
#                                                                     #
#  Example:  Using auxillary pair potentials with Tersoff's C-Si      #
#                                                                     #
#######################################################################
#                                                                     #
#    This file gives an example of how use auxillary pair potentials  #
#  with Tersoff's C-Si potential.  The energies and forces due to     #
#  these pair potentials are summed with the energy and forces due    #
#  to the Tersoff C-Si potential.  Pair potentials may be specified   #
#  between one or more pairs atom types - each such pair of atoms     #
#  will require its own potential.                                    #
#                                                                     #
#    The pair potential is entered as a polynomial of the form        #
#                                                                     #
#                       2         3                n+1                #
#          f(r) =   a  u   +  a  u   + ... + a    u                   #
#                    2         3              n+1                     #
#                                                                     #
#  where u = rc-r.  This polyomial form provides a simple and         #
#  compact specification while insuring that the potential            #
#  has a smooth cutoff of the r=rc.                                   #
#                                                                     #
#    The syntax of the pair potential command is                      #
#                                                                     #
#      POTENTIAL PAIR t1 t2 n rc a2 [a3 .. an+1]                      #
#                                                                     #
#    where                                                            #
#                                                                     #
#      t1,t2  are the atom types.                                     #
#      n      is the number of polymial terms.                        #
#      rc     is the pair potential cutoff in units of angtroms.      #
#      a2     is the first term, it is the coefficient of u^2.        #
#             It's units are (eunit/cm^2) where eunit is specified    #
#             by the EUNIT command.  The default value for eunit      #
#             is Kelvin.                                              #
#      a3     is the optional second term, the coefficient of u^3     #
#             It's units are (eunit/cm^3).                            #
#      an+1   is the optional n'th term, the coefficient of u^(n+1)   #
#             It's units are (eunit/cm^(n+1)).                        #
#                                                                     #
#######################################################################

#  Switch on Tersoff C-Si potential
potential set tersoff

#  Use old neighbor search algorithm to accomodate small atomic array
nsearch sort

#  Set up some useful constants for this potential
calc  C=1
calc  Si=2
calc  Al=3
calc  MassC=12.01
calc  MassSi=28.086
calc  MassAl=26.982
calc  A0=4.32
calc  DTIME=0.8e-15

#  Add pair potentials
calc NCOEFF=1
calc CUTOFF=3.0
calc COEFF=5e18
potential pair C  Al NCOEFF CUTOFF COEFF
potential pair Si Al NCOEFF CUTOFF COEFF
potential pair Al Al NCOEFF CUTOFF COEFF



#  Form diamond unit cell
box 2 2 2
particle  8
 C   0.125 0.125 0.125
 C   0.125 0.625 0.625
 C   0.625 0.125 0.625
 Al  0.625 0.625 0.125
 Si  0.375 0.375 0.375
 Si  0.375 0.875 0.875
 Si  0.875 0.375 0.875
 Al  0.875 0.875 0.375
dup  1  1 0 0
dup  1  0 1 0
dup  1  0 0 1

#  Scale to the unit cell length of diamond CSi
scale 4.32

#  Set particle masses
select type C
mass MassC
select type Si
mass MassSi
select type Al
mass MassAl

#  Set simulation step size
dtime DTIME

#  Instruct XMD to save energy every CMD step in file csi.e
esave 1 csipair.e

#  Initialize velocities to 1000K
itemp 1000

#  Set clamp off (adiabatic simulation)
clamp off

# Equilibrate particles for 500 step
write cor +csipair.cor
repeat 50
   cmd 10
   write cor +csipair.cor
end