#!/bin/sh
# tricking... the line after a these comments are interpreted as standard shell script \
    exec $ESPRESSO_SOURCE/Espresso $0 $*
#############################################################
#                                                           #
#  Lennard Jones Liquid                                     #
#                                                           #
#############################################################
#
# Copyright (C) 2010 The ESPResSo project
# Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany
#  
# This file is part of ESPResSo.
#  
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#  
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#  
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>. 
#  
#puts " "- just to put output to console
#set <var name> <value> - command to set a variable a value
    #a value can be constant or can be calculated one
    #a constant string value should be put with in " "
    #a number can be put with out " "
#a calculation can be done with in [ ]
    #Eg1: [analyze  mindist] : both are keywords..
    #Eg2: [expr 5*10]	     : expr used to do arithmetic
    #Eg3: [setmd time]	     : setmd is powerful command which can do lot of calculations
puts " "
puts "======================================================="
puts "=       dipole.tcl                                 ="
puts "======================================================="
puts " "

puts "Program Information: \n[code_info]\n"

#############################################################
#  Parameters                                               #
#############################################################

# System identification: 
set name  "dipole"
set ident "_test"

# 
set vmd_output "no"

# System parameters
#############################################################

# 10 000  Particles
set n_part  500
set density 0.1
set temperature 2.0


##############################################################
#SET UP THE THERMOSTAT
set gamma 1.0
thermostat langevin $temperature $gamma

##############################################################
#SET UP THE INTEGRATOR
setmd time_step 0.003

setmd skin      0.4
#############################################################
#SET UP PARTICLES
set volume [expr $n_part/$density]
set box_l [expr pow($volume, 1.0/3.0)]

setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1

#generate $n_part particles at random position
for {set i 0} { $i < $n_part } {incr i} {
    set posx [expr $box_l*[t_random]]
    set posy [expr $box_l*[t_random]]
    set posz [expr $box_l*[t_random]]

    set dipx [expr $box_l*[t_random]]
    set dipy [expr $box_l*[t_random]]
    set dipz [expr $box_l*[t_random]]
 
    part $i pos $posx $posy $posz dip $dipx $dipy $dipz dipm 1.0 type 0 mass 1.0
}

puts "Simulate System of $n_part particlesat density $density"
puts "Simulation box: $box_l"

#############################################################
# SET UP INTERACTION
set lj1_eps     1.0
set lj1_sig     1.0
set lj1_cut     1.1225
set lj1_shift   [calc_lj_shift $lj1_sig $lj1_cut]
set lj_offset   0.0

inter 0 0 lennard-jones $lj1_eps $lj1_sig $lj1_cut $lj1_shift $lj_offset

#############################################################
# MAGNETIC DIPOLE INTERACTION
set mdlc 2
set accuracy 1
set gap_size 1.1225

inter magnetic $mdlc $accuracy $gapsize[far_cutoff]
###############################################################
# WARM UP INTEGRATION
puts "\nWarmup:"

#after warmup, all particle pairs should have atleast this distance.
set min_dist     0.85

#compute minimum distance between two particles.
set act_min_dist [analyze mindist]

#set LJ capping
set cap 20
inter ljforcecap $cap

for {set i 0} { $act_min_dist < $min_dist } {incr i} {
    integrate 100

    #Warmup criterion
    set act_min_dist [analyze mindist]

    # Increase LJ cap
    set cap [expr $cap+10]
    inter ljforcecap $cap
}

# turn off capping
inter ljforcecap 0

###############################################################################
# MAIN INTEGRATION
puts "\n Main integration"

#open energy file
set dat_file [open "$name$ident-energy.dat" "w"]
puts $dat_file "\#md\ttmd\tE_tot\tE-kin\tE_LJ\tPressure\tE_mag"

# open VTF file
set vtf_file [open "$name$ident-conf.vtf" "w"]
writevsf $vtf_file

for {set i 0} {$i < 10 } {incr i } {
    integrate 10
    set tmd   [expr $i*0.003]
    set e_tot [analyze energy total]
    set e_kin [analyze energy kinetic]
    set e_lj [analyze energy nonbonded 0 0]
    set P_r [analyze pressure total]
    set e_mag [analyze energy magnetic]
  
    # output energies
    puts "\tstep=$i, time=[setmd time_step], e_tot=$e_tot, e_kin=$e_kin, e_lj=$e_lj, P_r=$P_r, e_mag=$e_mag"
    puts $dat_file "$i\t$tmd\t\t$e_tot\t\t$e_kin\t$e_lj\t$P_r\t$e_mag"

    # output configuration
    writevcf $vtf_file

    # store the last 20 configuarations
    #analyze push 20
}

close $dat_file
close $vtf_file

# terminate program
puts "\n\n Finished"