# we will use the following variable to calculate and print the average 
# number of ions found outside the sphere at each step

set avgNumIons 0
set Salto 0

wrapmode cell

##################################################

proc calcforces {step unique} {

  global Radius ForceConstant avgNumIons Salto max min delX
  
  if { $step > 0 && $step % 10000 == 0 } { 
	  
    set avgNumIons [expr $avgNumIons / 100.]
    print "Step $step, average number of ions outside the sphere: $avgNumIons"
	
    set avgNumIons 0
	  
    cleardrops 
	  
  }

  set ForceStep 0
  set IonsOutside 0
	  

	 while {[nextatom]} {

 
	  if { [getid] > 380 && [getid] < 38463 } { #se è acqua dimenticala
	      dropatom ;
	      incr Salto
	      continue
	    }


	    if { $step % 1 == 0 } {
	    
		    # vector between the ion and the sphere's center
		    
		    	set rvec [getcoord]
			foreach {x y z} $rvec { break } 		
			set absX [expr (abs($x))]


		    if { $absX < $min || $absX > $max} {
			#calcola di quanto è "fuori" lungo x rispetto  volume permesso in modo che poi gli applico in uno step una 				#forza sufficiente a farla rientrare
			#set delX [expr ($absX - $Radius)]
			#set deltaX [expr (abs($delX))]
			#set forceX [expr (-$ForceConstant * $deltaX)]
			if { $x < $min} {
			set delX [expr ($absX - $min)]
			set deltaX [expr (abs($delX))]
			set forceX [expr ($ForceConstant * $deltaX)]
			addforce "$forceX 0 0"
			}
			if {$x > $max} {
			set delX [expr ($absX - $max)]
			set deltaX [expr (abs($delX))]
			set forceX [expr (-$ForceConstant * $deltaX)]
			addforce "$forceX 0 0"
			}
			set ForceStep [expr ($ForceStep + abs($forceX))]
			print "Step $step 	Force $forceX, Tot $ForceStep"
			incr avgNumIons
			incr IonsOutside
			#} else {
			#  dropatom ;# this ion is already inside the sphere
		    }
	      }
	   }
  
  # at the end of cycle over atoms print info every n steps

  if { $step % 1000 == 0 } {
	 print "Step $step	Force $ForceStep   IonsOutside $IonsOutside"  
         print "atomi saltati $Salto"
  }

}