######################### BEGIN USER-DEFINED THINGIES ##################################
# From Ballenegger: T = 300 K;  mu = 2.45 D ; sigma = 0.3668 nm ; u = 2.1747 kJ/mol ; 
# Reduced quantities are rho* = rho sigma^3 ; T* = kT / u ; mu* = sqrt[mu^2/(sigma^3 u)] 
# TD state is rho*=0.82 ; T*=1.15 ; mu*=1.82
########################################################################################


########################################################################################
########################## Constants Used in the Simulation  ###########################
########################################################################################
#Length of the Dipole=
set red_diplen 0.3   
#Number of equilibration steps=                           
set eqsteps 5000
#Number of sample steps=
set sampsteps 16000000
#Total number of Dipoles=
set npart 125
#Density=
set red_rho 0.82
#Temperature=
set red_T 1.15
#Dipolar Moment=
set red_mu 1.82 
#Cutoff distance for the Lennard-Jones potential=
set red_cutoff [expr 0.8 / 0.3668 ] 
#Time Step=
setmd time_step 0.01
#Skin for MD=
setmd skin 0.3
#PI constant=
set PI 3.14159265358979323844



#######################################################################################
########################## Declaration of Variables ###################################
#######################################################################################
set mol_id 0
set time_step [setmd time_step]

#Set Particle Type 
set typecounter 0
set LJ $typecounter ; incr typecounter
set POS $typecounter ; incr typecounter
set NEG $typecounter ; incr typecounter

#Set Size of the Simulation Box
set red_vol [expr $npart / $red_rho ]
set red_boxl [expr pow($red_vol,1./3.)]
setmd box_l $red_boxl $red_boxl $red_boxl

#Set Interaction Type 
set intcounter 0
set intLJ $intcounter ; incr intcounter
set intCONSTR_QQ $intcounter ; incr intcounter
set constr_qq $red_diplen
#inter $intCONSTR_QQ rigid_bond $constr_qq    0.01 0.01
inter $LJ $LJ lennard-jones 1.0 1.0 $red_cutoff 0. 0.


#Set Charge of Positive and Negative Particles
set chargep [expr $red_mu / $red_diplen ]
set chargem [expr -$chargep ]

#Set Thermostat
thermostat off
thermostat langevin $red_T [expr 1./(25.*$time_step)]


#Function to calculate the Total Dipole Moment of the System 
proc dipolemoment {} {
        global chargep
        set totdip {0 0 0} 
        set dip {0 0 0}
	for { set i 0} { $i < [setmd n_part] } { set i [expr $i+3] } { 
                set dv [bond_vec_min [part [expr $i+1]  pr pos] [part [expr $i+2] pr pos]] 
	        set dq $chargep
	        set dip [vecscale $dq $dv]
                set totdip [vecadd $totdip $dip]
        }   
        return  $totdip
}

#Function to locate the Particles in the Simulation Box 
proc create_mol { posx posy posz orx ory orz } {
        global mol_id
        global red_diplen
        global LJ
        global POS
        global NEG
        global chargep
        global chargem
        global intCONSTR_QQ
        global particle_list

        set d [expr $red_diplen/2.]
        set len [expr sqrt($orx*$orx+ $ory*$ory + $orz*$orz)]
        set orx [expr $orx/$len] ; set ory [expr $ory/$len] ; set orz [expr $orz/$len]

        #Central Particle
        set LJid [setmd n_part]
        part $LJid pos $posx $posy $posz type $LJ q 0 virtual 0 mol $mol_id mass 1 

        #Positive Particle
        set POSid [setmd n_part]
        part $POSid pos [expr $posx + $d*$orx ] [expr $posy + $d*$ory] [expr $posz + $d*$orz ]  type $POS q $chargep virtual 1 mol $mol_id mass 1

        #Negative Particle
        set NEGid [setmd n_part]
        part $NEGid pos [expr $posx - $d*$orx ] [expr $posy - $d*$ory] [expr $posz - $d*$orz ]  type $NEG q $chargem virtual 1  mol $mol_id mass 1

#        part $POSid bond $intCONSTR_QQ $NEGid
        set pair [list $mol_id $LJid $POSid $NEGid ] 
#        set pair [list $mol_id $LJid  ] 
        lappend particle_list $pair 
        incr mol_id
}

puts "Placing Particles..."

#Loop to Create Molecules Uniformly Spaced
set x 0
set y 0
set z 0
set nmax [expr int($red_boxl)+1 ]
puts "max number of cells $nmax"
set indic 0
for { set i 1} { $i < $nmax } { incr i} {
     set x [expr (1.*$i-1.0) + 0.5 ]
     set y 0.5
     for { set j 1} { $j < $nmax } { incr j} {
            set y [expr (1.*$j-1.0) + 0.5 ]
            set z 0.5
            for { set k 1} { $k < $nmax } { incr k} {
                  set z [expr (1.*$k-1.0) + 0.5 ]

                  if { $indic  < 2  } {
                  incr indic
                     create_mol $x $y $z  [t_random] [t_random] [t_random]
                  }
            }
     }
}
 

puts "number of particles $indic"
eval analyze set $particle_list
puts "[ analyze set topo_part_sync $particle_list ]"
set particles [setmd n_part]
puts "Number of particles $particles "
puts ""
puts "Finished Placing Particles"


#Where to Write to
set dir "/auto.anoa/home/pojeda/DIPOLAR_LIQUIDS/"
# Tagging of Output Files
set tag [format "$dir/VS_W-Q-%.2f-mu-%.2f-rho-%.2f-tmp-%.2f-tstp-%.2f-lenght-%.2f" $chargep $red_mu $red_rho $red_T $time_step $red_diplen]



#Write the Initial Configuration to a PDB File
writepdb "inicial.pdb"


#Tuning P3M Algorithm
puts "verlet reuse [setmd verlet_reuse]  max_range [setmd max_range]" 
puts "Tuning P3M..."
puts "[inter coulomb [expr 1.0/$red_T ] p3m tunev2 accuracy 1e-2 r_cut $red_cutoff  mesh 32 ]"
puts "verlet reuse [setmd verlet_reuse]  max_range [setmd max_range]" 


for { set i 1 } { $i < 100 }  { incr i } {
       integrate 1
       puts "$i [analyze energy coulomb]"       
}



exit

set MM 0
#Equilibration Phase
puts "Equilibration phase"
for { set i 1 } { $i < 100 }  { incr i } {
       integrate 1
       set dipmom [dipolemoment]
       set Mx [lindex $dipmom 0] ; set My [lindex $dipmom 1] ; set Mz [lindex $dipmom 2] ;
       set instantMM [expr $Mx*$Mx+$My*$My+$Mz*$Mz]
       set MM [expr $MM + $instantMM]  ;
       set avgMM [ expr $MM/$i] ;
       set eps [expr 1+3*(4.0*$PI/$red_T)*($red_rho/9.)/$npart * $avgMM] ;
#       puts "$i $eps"
       puts "$i [analyze energy coulomb]"
       writepdb "file.pdb" 
}


exit