# 3D Gay Berne on Surface simulation
units real
atom_style ellipsoid
dimension 3
boundary p p p
atom_modify id yes
atom_modify map array
atom_modify sort 100 25
lattice fcc 5.0
region box block 0.0 5.0 0.0 5.0 0.0 5.0
create_box 1 box
#create_atoms 1 region substrate
create_atoms 1 single 2 2 2
#create_atoms 1 single 3 3 3
compute_modify extra 0
set type 1 mass 16.04
set type 1 shape 1 1 1
velocity all create 300.0 4928459 rot yes dist gaussian
set group all quat/random 18238
compute rot all temp/asphere
compute_modify extra 0
group spheroid type 1
variable dof equal count(spheroid)+2
compute_modify rot extra ${dof}
neigh_modify delay 0 every 1
timestep 0.0002
variable z_atom atom z
variable x_atom atom x
variable y_atom atom y
# Get Constants of the Energy That Don't Change Upon Taking Directional Derivatives
variable energy_x atom (1.47025e-06*v_x_atom^4+-0.000598126*v_x_atom^2+0.0156679+0.00012389*v_x_atom^3+-0.0128309*v_x_atom)*sin(1.04344*v_x_atom)+-1098.47+sin(-10.7193*v_x_atom+0.464616)*0.000818385+0.24923*cos(0.0319159*v_x_atom+0.106606)
variable energy_y atom (1.47025e-06*v_y_atom^4+-0.000598126*v_y_atom^2+0.0156679+0.00012389*v_y_atom^3+-0.0128309*v_y_atom)*sin(1.04344*v_y_atom)+-1098.47+sin(-10.7193*v_y_atom+0.464616)*0.000818385+0.24923*cos(0.0319159*v_y_atom+0.106606)
variable energy_z atom 6602.67*exp(-2.41063*v_z_atom)-176.047*v_z_atom^(-3)-1829.08*v_z_atom^(-6)+2866.59*v_z_atom^(-8)
variable energy_total atom v_energy_x*v_energy_y*v_energy_z
#Trans Forces
variable fx_2 atom v_energy_y*v_energy_z*(-5.881*10^(-6)*v_x_atom^3+1.19624*10^(-3)*v_x_atom-3.7167*10^(-4)*v_x_atom^2+.0128309)*sin(1.04344*v_x_atom)+(-1.47025*10^(-6)*v_x_atom^(4)+.00059812*v_x_atom^2-.0156679-.00012399*v_x_atom^3+.0128309*v_x_atom)*1.04344*cos(1.04344*v_x_atom)+8.772139*10^(-3)*cos(-10.713*v_x_atom+.464616)+7.954399*10^(-3)*sin(.0319159*v_x_atom+.106606)
variable fy_2 atom v_energy_x*v_energy_z*(-5.881*10^(-6)*v_y_atom^3+1.19624*10^(-3)*v_y_atom-3.7167*10^(-4)*v_y_atom^2+.0128309)*sin(1.04344*v_y_atom)+(-1.47025*10^(-6)*v_y_atom^(4)+.00059812*v_y_atom^2-.0156679-.00012399*v_y_atom^3+.0128309*v_y_atom)*1.04344*cos(1.04344*v_y_atom)+8.772139*10^(-3)*cos(-10.7193*v_y_atom+.464616)+7.954399*10^(-3)*sin(.0319159*v_y_atom+.106606)
variable fz_2 atom v_energy_x*v_energy_y*15917.0*exp(-2.41063*v_z_atom)-528.14*v_z_atom^(-4)+10974*v_z_atom^(-7)+22932.7*v_z_atom^(-9)
compute q all property/atom quatw quati quatj quatk
fix 1 all nve
#fix 1 all nvt temp 300.0 300.0 1.0
dump 1 all custom 100000 dump.ellipse.gayberne &
id type x y z c_q[1] c_q[2] c_q[3] c_q[4]
dump 2 all image 100000 image.*.jpg type type adiam 1.2
dump_modify 2 pad 5
dump 4 all custom 100000 dump.force.test &
id v_fx_2 v_fy_2 v_fz_2 v_energy_x v_energy_y v_energy_z
dump 5 all custom 100000 dump.x.test &
id x
thermo 10000
thermo_style custom step temp ke pe etotal press vol
run 50000000
--