# Script to update simulation parameters (new target coordinates, spring coefficients etc.) and save them replacing the old ones)

# ______________ User parameters section  ==>
# These are default values.
# Any/all of these parameters can be overwritten by the same parameters defined in configuration file.

# ==>	Working folder
set workingFolder {/simulations/test_steered_simulation}

# ==>	Configuration file name
set configFileName update_parameters_conf.txt




# _________________________ Computational section
set start_time [clock seconds]

set oldFolder [pwd]
cd $workingFolder

# ___________________ reading parametrization parameters from configuration file
# Open configuration file
set fileLook [open $configFileName r]
# Read configuration file and recognize data
if {[file exists $configFileName] == 0} {
	# Configuration file does not exist. Stopping the minimization...
	puts "Configuration file does not exist. Stopping the minimization..."
} else {
	# Configuration file does exist
	puts "Reading configuration file $configFileName"
	while {[gets $fileLook line] >= 0} {
		set lineRecognize [lindex $line 0]
		#puts "'$line'	'$lineRecognize'"
		switch -regexp $lineRecognize {
			"^$" {}
			"^#" {}
			"^workingFolder$" {
				puts "$line"
				set workingFolder [lindex $line 1]
				cd $workingFolder
			}
			"^configFileName$" {
				# Configuration file name
				puts "$line"
				set configFileName [lindex $line 1]
			}
			"^statusFileName$" {
				# Status file name
				puts "$line"
				set statusFileName [lindex $line 1]
			}
			"^structure$" {
				# Protein structure file
				puts "$line"
				set structure [lindex $line 1]
			}
			"^coordinates$" {
				# Protein coordinates
				puts "$line"
				set coordinates [lindex $line 1]
			}
			"^coordinatesFileType$" {
				# Coordinates File type
				puts "$line"
				set coordinatesFileType [lindex $line 1]
			}
			"^referenceCoordinates$" {
				# Reference coordinates - used for various calculations as a steady reference point
				puts "$line"
				set referenceCoordinates [lindex $line 1]
			}
			"^referenceCoordinatesFileType$" {
				# Reference Coordinates File type
				puts "$line"
				set referenceCoordinatesFileType [lindex $line 1]
			}
			"^constantsFileName$" {
				# Output spring constant coefficients file name, pdb type only
				puts "$line"
				set constantsFileName [lindex $line 1]
			}
			"^constantsColumn$" {
				# Input spring constant coefficients column
				puts "$line"
				set constantsColumn [string tolower [lindex $line 1]]
			}
			"^outputCoordinatesFileName$" {
				# Output coordinates file name
				puts "$line"
				set outputCoordinatesFileName [lindex $line 1]
			}
			"^outputConstantsFileName$" {
				# Output spring constant coefficients file name
				puts "$line"
				set outputConstantsFileName [lindex $line 1]
			}
			"^timeAveraging$" {
				#  Option to average among time frames 1 - yes, 0 - no 
				puts "$line"
				set timeAveraging [lindex $line 1]
			}
			"^realFrame$" {
				#  Option to choose the real frame closest to an average (by rmsd) as an average frame   1 - yes, 0 - no
				puts "$line"
				set realFrame [lindex $line 1]
			}
			"^saveTimeAveragedStructure$" {
				#  Option to Save time-averaged structure
				puts "$line"
				set saveTimeAveragedStructure [lindex $line 1]
			}
			"^selectionText$" {
				#  Selection for atoms which should be parametrized
				puts "$line"
				set selectionText [lindex $line 1]
			}
			"^firstFrame$" {
				#   start frame for time averaging. If < 0, then this frame is relative to the last frame
				puts "$line"
				set firstFrame [lindex $line 1]
			}
			"^lastFrame$" {
				#   Last frame for time averaging. If -1, then the very last frame
				puts "$line"
				set lastFrame [lindex $line 1]
			}
			"^preserveUnassignedConstants$" {
				# If <1>, Read spring constants from the existing file, calculate new constants, assign them and leave all old constants, which were not reassigned, in place. Useful when there are some additional harmonic restraints in the system, not related to parametryzation. If <0> - set all the coefficients except newly calculated, to 0.
				puts "$line"
				set preserveUnassignedConstants [lindex $line 1]
			}
			"^preserveUnassignedTargets$" {
				# If <1>, Read Target coordinates from the existing file, calculate new target coordinates, assign them and leave all old coordinates, which were not reassigned, in place. Useful when there are some additional harmonic restraints in the system, not related to parametryzation. If <0> - set all the coordinates except newly calculated, to 0. 
				puts "$line"
				set preserveUnassignedTargets [lindex $line 1]
			}
			"^independentConstants$" {
				# <1> - Independent initial spring constants for every atom. Then they scale according to constantsScaleMode. <0> - equal initial spring constants for all the restrained atoms. Then they scale according to constantsScaleMode. 
				puts "$line"
				set independentConstants [lindex $line 1]
			}
			"^targetCalculationMode$" {
				# Mode of the target coordinates calculation
				# <0> - "fixed target" mode: no target update
				# <1> - "cylindrical projection target": target points are mapped onto cyllinder of predefined (in u"pdate_parameters.vmd" script) radius. They have the same angular coordinates as in the input coordinates file
				# <2> - "constant force outward": all the restrained atoms get the same amount of force in the direction from the central (z) axis through the atom.
				# <3> - "Coriolis force": all the restrained atoms get the force propotional to the square of the distance to the central axis. Force will be in the direction from the central (z) axis through the atom
				# <4> - "constant shift": target coordinate of the atom shifts every time in radial direction by the fixed distance relative to the current coordinate of the atom
				# <5> - "cylindrical projection (selection 1) and axis (selection 2) target" : target points are mapped onto cyllinder of predefined (in the "update_parameters.vmd" script) radius. They have the same angular coordinates as in the input coordinates file
				# <6> - "constant force outward (selection 1) and axis (selection 2)" : all the restrained atoms get the same amount of force in the direction from the central (z) axis through the atom.
				# <7> - "Coriolis force outward (selection 1) and to the axis (selection 2)" : all the restrained atoms get the force propotional to the square of the distance to the central axis. Force will be in the direction from the central (z) axis through the atom
				# <8> - "constant shift outward (selection 1) and to the axis (selection 2)" : target coordinate of the atom shifts every time in radial direction by the fixed distance relative to the current coordinate of the atom
				# <9> - "restrained center of mass of every subselection" : atoms of every subselection provided as a sublist in the "selection" variable (for one selection use one sublist like {{resname ALA}} ) will be harmonically restrained to tranlsated location such that COM of every subselection will be restrained to it's location in the reference file
				# <10> - "restrained center of mass of every fragment in selection" : atoms of every fragment in the "selection" will be harmonically restrained to tranlsated location such that COM of every fragment will be restrained to it's location in the reference file

				puts "$line"
				set targetCalculationMode [lindex $line 1]
			}
			"^constantsScaleMode$" {
				# Mode of scaling of the spring constants: 
				# <0> - 'fixed constant'. 
				# <1> - 'exponential increase'. 
				# <2> - 'constant force' mode (constants adjusted such that shift by 1 A gives initially the same amount of energy. 
				# <3> - 'RMSD gradient' mode: spring constant is increased $multiplierConstant times, if relative decrease in RMSD during the last simulation step was smaller, than required value $RmsdGradient (particular case can be RmsdGradient = 0, which maintains the simulation at constant RMSD from symmstric structure); If RMSD decrease was larger, than required $RmsdGradient amount, then spring constant is decreased $multiplierConstant times.
				# <4> - "oscillating constant" : spring constant amount will periodically increase and decrease exponentially (linear in log scale) between the $oscillatingConstantMin and $oscillatingConstantMax with the rate of $multiplierConstant/step and $oscillatingConstantPeriod (once in a oscillating period it will go up and down). The first hal-cycle will be decrease of the constant, so initial constants should be at their maximal values. After complete reductin cycle they should not be less than 0.001
				# <5> - "linear increase". Coefficient=initialConstant+multiplierConstant*timestep. For "smooth" system should lead to the constant force mode.
				# <6> - "parabolic increase". Coefficient=initialConstant+(multiplierConstant*timestep)^2. For "smooth" system should lead to the constant energy per A shift mode.
				# <7> - "square-hyperbolic increase" Coefficient=initialConstant/(1-multiplierConstant*timestep)^2. For "smooth" system should lead to the linear RMSD decrease mode. multiplierConstant*timestep should be < 1. multiplierConstant = (1-initialConstant/maxconstant)/timesteps
				puts "$line"
				set constantsScaleMode [lindex $line 1]
			}
			"^initialConstant$" {
				# Initial value of harmonic spring constant
				puts "$line"
				set initialConstant [lindex $line 1]
			}
			"^multiplierConstant$" {
				# This is the multiplier for exponential growth of the spring constant
				puts "$line"
				set multiplierConstant [lindex $line 1]
			}
			"^energySlopeConstant$" {
				# Amount of energy per one angstrom for "constant force" coefficient mode
				puts "$line"
				set energySlopeConstant [lindex $line 1]
			}
			"^RmsdGradient$" {
				# Target ratio of new RMSD to previous RMSD (for constantsScaleMode	3 "RmsdGradient" method of constant scale mode)
				puts "$line"
				set RmsdGradient [lindex $line 1]
			}
			"^oscillatingConstantMin$" {
				# Minimal value of spring constant for "oscillating constant" spring change mode
				puts "$line"
				set oscillatingConstantMin [lindex $line 1]
			}
			"^oscillatingConstantMax$" {
				# Maximal value of spring constant for "oscillating constant" spring change mode
				puts "$line"
				set oscillatingConstantMax [lindex $line 1]
			}
			"^oscillatingConstantPeriod$" {
				# Oscillation period (in steps) of spring constant for "oscillating constant" spring change mode.
				puts "$line"
				set oscillatingConstantPeriod [lindex $line 1]
			}
			"^targetRmsd$" {
				# Maximal acceptable target value of RMSD from the average structure. If this value (or lower) was reached, simulation stops
				puts "$line"
				set targetRmsd [lindex $line 1]
			}
			"^parametrizationStepsMax$" {
				# Maximal possible number of parametrization steps, before it stops (if was not stopped earlier by other criterion).
				puts "$line"
				set parametrizationStepsMax [lindex $line 1]
			}
			"^logOutputFiles$" {
				# Option to keep log (numbered copies) of all output files - parametrized structure, constants file, status file
				puts "$line"
				set logOutputFiles [lindex $line 1]
			}
			"^logOutputFormat$" {
				# Format to keep log of structure-format output files (target coordinates and spring constants)
				puts "$line"
				set logOutputFormat [lindex $line 1]
			}
			"^cylinderRadius$" {
				# Target Cylinder radius for "cylindrical projection target" calculation mode
				puts "$line"
				set cylinderRadius [lindex $line 1]
			}
			"^constantShift$" {
				# Shift per parametrizatin step for "constant shift" mode: target coordinate of the atom shifts every time in radial direction by the fixed distance relative to the current coordinate of the atom
				puts "$line"
				set constantShift [lindex $line 1]
			}
			"^twist$" {
				# Twist (in degrees) for the direction of force for different restrain modes
				puts "$line"
				set twist [lindex $line 1]
			}
			"^selection2Text$" {
				# Second Selection for atoms which should be dynamically restrained - i.e. for which all these calculations are performed. It is used in complex restraints
				puts "$line"
				set selection2Text [lindex $line 1]
			}
			"^twist2$" {
				# Twist (in degrees) for the direction of force for different restrain modes. Positive direction is counterclockwise looking from the top of z axis. Parameter for the selection2 for complex restraints
				puts "$line"
				set twist2 [lindex $line 1]
			}
			default {
				puts "Unrecognized parameter   '$lineRecognize'"
			}
		}
	}
	close $fileLook
	
	# ________ Autocalculate the missing constants
	if {[llength $multiplierConstant]==0} {;	# Autocalculate multiplierConstant
		# Mode of scaling of the spring constants : 
		case $constantsScaleMode {
			"0" {
				# <0> - "fixed constant"
				set multiplierConstant 1.0
			}
			"1" {
				# <1> - "exponential increase" . 
# 				set multiplierConstant 9999.999/($initialConstant*$parametrizationStepsMax)
				if {([catch {set $oscillatingConstantMax}]==0)&&([llength $oscillatingConstantMax]>0)} {
					set multiplierConstant [expr {pow(1.0*$oscillatingConstantMax/$initialConstant,1.0/$parametrizationStepsMax)}]
				} else {
					set multiplierConstant [expr {pow(10.0/$initialConstant,1.0/$parametrizationStepsMax)}]
				}
			}
			"2" {
				# <2> - "constant force" mode (constants adjusted such that shift by 1 A gives initially the same amount of energy. 
				set multiplierConstant 2
			}
			"3" {
				# <3> - "RMSD gradient" mode : spring constant is increased $multiplierConstant times, if relative decrease in RMSD during the last simulation step was smaller, than required value $RmsdGradient (particular case can be RmsdGradient = 0, which maintains the simulation at constant RMSD from symmstric structure); If RMSD decrease was larger, than required $RmsdGradient amount, then spring constant is decreased $multiplierConstant times.
				set multiplierConstant 2
			}
			"4" {
				# <4> - "oscillating constant" : spring constant amount will periodically increase and decrease exponentially (linear in log scale) between the $oscillatingConstantMin and $oscillatingConstantMax with the rate of $multiplierConstant/step and $oscillatingConstantPeriod (once in a oscillating period it will go up and down). The first hal-cycle will be decrease of the constant, so initial constants should be at their maximal values. After complete reductin cycle they should not be less than 0.001
				if {([catch {set $oscillatingConstantMax}]==0)&&([llength $oscillatingConstantMax]>0)} {
					set multiplierConstant [expr {3.0*$oscillatingConstantMax/($initialConstant*$parametrizationStepsMax)}]
				} else {
					set multiplierConstant [expr {30.0/($initialConstant*$parametrizationStepsMax)}]
				}
			}
			"5" {
				# <5> - "linear increase" Coefficient=initialConstant+multiplierConstant*timestep. For "smooth" system should lead to the constant force mode.
				if {([catch {set $oscillatingConstantMax}]==0)&&([llength $oscillatingConstantMax]>0)} {
					set multiplierConstant [expr {double($oscillatingConstantMax-$initialConstant)/$parametrizationStepsMax}]
				} else {
					set multiplierConstant [expr {(10.0-$initialConstant)/$parametrizationStepsMax}]
				}
			}
			"6" {
				# <6> - "parabolic increase" Coefficient=initialConstant+(multiplierConstant*timestep)^2. For "smooth" system should lead to the constant energy per A shift mode.
				if {([catch {set $oscillatingConstantMax}]==0)&&([llength $oscillatingConstantMax]>0)} {
					set multiplierConstant [expr {($oscillatingConstantMax-$initialConstant)/pow($parametrizationStepsMax,2)}]
				} else {
					set multiplierConstant [expr {(10-$initialConstant)/pow($parametrizationStepsMax,2)}]
				}
			}
			"7" {
				# <7> - "square-hyperbolic increase" Coefficient=initialConstant/(1-multiplierConstant*timestep)^2. For "smooth" system should lead to the linear RMSD decrease mode. multiplierConstant*timestep should be < 1. multiplierConstant = (1-initialConstant/maxconstant)/timesteps
					set multiplierConstant [expr {(1.0-$initialConstant/9999.999)/$parametrizationStepsMax}]
			}
			default {
				set multiplierConstant 1.0
			}
		};	# Mode of scaling of the spring constants 
		puts "Autocalculate multiplierConstant: $multiplierConstant"
	};	# Autocalculate multiplierConstant
}


set statusLogFileName [join "$statusFileName _log.txt" {}]

# _________________ reading parametrization status from the status file
# Check if status file exists
if {[file exists $statusFileName] == 0} {
	# Status file does not exist. Initialize status parameters
	puts "Status file does not exist. Initializing status parameters."
	set parametrizationCyclesCompleted 0
	set previousConstant $initialConstant
	
	# Make a copy for the log, if necessary
	if {$logOutputFiles == 1} {
		set fileLook [open $statusLogFileName w]
		puts $fileLook "# Autogenerated minimization status file log"
		puts $fileLook "parametrizationCyclesCompleted\tpreviousRmsd\tpreviousConstant"
		close $fileLook
	} 
} else {
	# status file exists. Read it.
	# This is not the first minimization cycle
	# open configuration file
	set fileLook [open $statusFileName r]
	# read configuration file and recognize data
	while {[gets $fileLook line] >= 0} {
		set lineRecognize [lindex $line 0]
		switch -regexp $lineRecognize {
			"^parametrizationCyclesCompleted$" {
				# Number of minimization cycles completed up to the moment
				puts "$line"
				set parametrizationCyclesCompleted [lindex $line 1]
				incr parametrizationCyclesCompleted
			}
			"^previousRmsd$" {
				# Average RMSD between the monomers and parametric targed, obtained in the last step
				puts "$line"
				set previousRmsd [lindex $line 1]
			}
			"^previousConstant$" {
				# Spring constant value attained in the previous step
				puts "$line"
				set previousConstant [lindex $line 1]
			}
			default {
				puts "Unrecognized parameter   '$lineRecognize'"
			}
		}
	}
	close $fileLook
}


# ___________________________ parametrization
if {($timeAveraging == 1) || ($targetCalculationMode!=0) || ($constantsScaleMode!=0) || ($parametrizationCyclesCompleted==0)} {
	# There will be some time-averaging or parametrization
	# Read the structure
	if {[catch {
		mol new $structure type psf
		mol addfile $coordinates type $coordinatesFileType waitfor all
		set coordinatesMolID [molinfo top get id]
	}]} {error "Can't load the structure"}
	puts "Structure loaded"

	# convert selection phrase into atomselection
	# If there will be 2 selectionz used for target calculation?
	case $targetCalculationMode {
		{0 1 2 3 4 10} {
			# There will be 1 selection (selection 2 will be ignored)
			set selection [atomselect top $selectionText]
		}
		{5 6 7 8} {
			# There will be 2 selections. The same Spring constants will be calculated for joined selections 1 and 2
			set selection1 [atomselect top $selectionText]
			set selection2 [atomselect top $selection2Text]
			set selection [atomselect top "($selectionText) or ($selection2Text)"]
		}
		{9} {
			# Every sublist of the selection text will be processed into one selection with "or" as a joint
			set selection [atomselect top [join $selectionText { or }]]
		}
		default {
			error "Unrecognized target calculation mode. Exiting."
		}
	};	# If there will be 2 selections used?

	# repositioned atoms will be marked by positive betas, all other - by zero beta
	set all [atomselect top all]
	$all set beta 0
	$selection set beta 1
	
	
	# List of indexes in the selection
	set selectionIndexes [$selection get index]
	
	# Define numbers of the first and last frames
	set framesNumber [molinfo top get numframes]
	if {$lastFrame == -1 || $lastFrame > [expr {$framesNumber-1}]} {set lastFrame [expr {$framesNumber-1}]}
	if {$firstFrame < 0} {
		# First frame is several steps back relative to the last frame
		set firstFrame [expr {$lastFrame + $firstFrame}]
		if {$firstFrame < 0} {set firstFrame 0}
	}
	
	# Duplicate the last frame - this new frame will get the calculated average values and will be used for comparison
	animate dup frame $lastFrame [molinfo top]
	set framesNumber [molinfo top get numframes]
	set averageFrame [expr {$framesNumber - 1}]
	$all frame $averageFrame
	$all update
}


# _______ make averaging among timesteps
if {$timeAveraging == 1} {
	puts "calculating average among timesteps..."
	# Go through the frames and average coordinates for selection
	
	# Update the frame
	$selection frame $firstFrame
	$selection update
	
	puts "frame $firstFrame of the range from $firstFrame to $lastFrame..."
	set averageX [$selection get x]
	set averageY [$selection get y]
	set averageZ [$selection get z]
	for {set frame [expr {$firstFrame+1}]} {$frame <= $lastFrame} {incr frame} {
		puts "frame $frame of the range from $firstFrame to $lastFrame..."
		# Update the frame
		$selection frame $frame
		$selection update

		set averageX [vecadd $averageX [$selection get x]]
		set averageY [vecadd $averageY [$selection get y]]
		set averageZ [vecadd $averageZ [$selection get z]]
	}
	set scalingFactor [expr {1.0/($lastFrame-$firstFrame+1)}]
	set averageX [vecscale $scalingFactor $averageX]
	set averageY [vecscale $scalingFactor $averageY]
	set averageZ [vecscale $scalingFactor $averageZ]
	
	
	if {$realFrame == 1} {
		# choose the real frame closest to an average (by rmsd) as an average frame or/and calculate the standard deviation of all the frames from the average
		puts "comparing all the timesteps with an average..."
		
		# Set averaged values to the new average frame
		$selection frame $averageFrame
		$selection update
		$selection set x $averageX
		$selection set y $averageY
		$selection set z $averageZ

		
		# Compare all other frames with the first/averaged frame
        set averageFrameSelection [atomselect top "index $selectionIndexes" frame $averageFrame]
		set rmsdMin -1
		for {set frame [expr {$firstFrame}]} {$frame <= $lastFrame} {incr frame} {
			puts "2nd pass: frame $frame of the range from $firstFrame to $lastFrame..."
			# Update the frame
			$selection frame $frame
			$selection update
			
			# Calculate an average rmsd for this whole frame
			if {$realFrame == 1} {
		        set rmsd [measure rmsd $selection $averageFrameSelection]
		        if {$rmsdMin < 0 || $rmsdMin > $rmsd} {
			        # the first comparison or smaller rmsd
			        set rmsdMin $rmsd
			        set rmsdMinFrame $frame
		        }
			}
		}
		
		if {$realFrame == 1} {
			# Record coordinates of the frame closest to average as the average coordinates
			$selection frame $rmsdMinFrame
			$selection update
			set averageX [$selection get x]
			set averageY [$selection get y]
			set averageZ [$selection get z]
		}
		
	}
	
	# Assign averaged values to the last frame
	$selection frame $averageFrame
	$selection update
	$selection set x $averageX
	$selection set y $averageY
	$selection set z $averageZ
	$all frame $averageFrame
	$all update
	
	if {$saveTimeAveragedStructure ==1} {
		# Save time-averaged structure
		$all writepdb [molinfo top get name]_time-avrg.pdb
		# Make a copy for the log, if necessary
		if {$logOutputFiles == 1} {
			# file copy -force $outputConstantsFileName $outputConstantsFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted]
			switch $logOutputFormat {
				"pdb" {set logOutputFormatExtension "pdb"}
				"dcd" {set logOutputFormatExtension "dcd"}
				"namdbin" {set logOutputFormatExtension "coor"}
				default {set logOutputFormatExtension "pdb"}
			}
					
			animate write $logOutputFormat [molinfo top get name]_time-avrg[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted].$logOutputFormatExtension beg $averageFrame end $averageFrame
		}
	}
} else {
}






# ____________ Calculating the new spring constants and target coordinates
if {($targetCalculationMode!=0) || ($constantsScaleMode!=0) || ($parametrizationCyclesCompleted==0)} {
	# Remember the current target coordinates - to restore them after the possible constants assignement into x,y,z columns
	set targetXYZ [$all get {x y z}]

	
	# _______________________ calculation of spring constant coefficients
	puts "Calculating and recording spring constants..."
	
	# Calculate RMSD between the all/last frame and the target !!!!!!!!!!!!!!!!!!!!!! Not implemented yet
    set averageFrameSelection [atomselect top "index $selectionIndexes" frame $averageFrame]
	if {$timeAveraging == 1} {
		set currentRmsd 0
		for {set frame [expr {$firstFrame}]} {$frame <= $lastFrame} {incr frame} {
			# Update the frame
			$selection frame $frame
			$selection update
			
			# Calculate an average rmsd for this whole frame
	        set rmsd [measure rmsd $selection $averageFrameSelection]
	        set currentRmsd [expr {$currentRmsd+pow($rmsd,2)}]
        }
        set currentRmsd [expr {pow(($currentRmsd/($lastFrame-$firstFrame+1)),0.5)}]
	} else {
		$selection frame $lastFrame
		$selection update
		set currentRmsd [measure rmsd $selection $averageFrameSelection]
	}
	
	# Read spring constants from file, if necessary
	if {$preserveUnassignedConstants ==1 || $independentConstants == 1} {
		# Read spring constants from the existing file, calculate new constants, assign them and leave all old constants, which were not reassigned, in place
		# Useful when there are some additional harmonic restraints in the system, not related to parametryzation
		if {[catch {
			# Try to read old spring constants
			mol new $structure type psf
			mol addfile $outputConstantsFileName type pdb waitfor all
		}] == 0} {
			# Old spring constants were read
			$all set $constantsColumn [[atomselect top all] get $constantsColumn]
			mol delete top
		} else {
			# There was not constants file. Initialize one
			# Zero the whole column for the future spring constant values. Set predefined constant for restrained selection
			if {[molinfo top get numframes]==0} {
				# unsuccessfully loaded PSF
				mol delete top
			}
			$all set $constantsColumn 0
			$averageFrameSelection set $constantsColumn $initialConstant
		}
		set previousConstants [$averageFrameSelection get $constantsColumn]
	} else {
		# All old sonstants are set to 0
		$all set $constantsColumn 0
		$averageFrameSelection set $constantsColumn $initialConstant
	}		
	
	
	# ___ Calculate new spring constants
	# Mode of scaling of the spring constants: <0> - fixed constant. <1> - exponential increase. <2> - "constant force" mode (constants adjusted such that shift by 1 A gives initially the same amount of energy
	switch $constantsScaleMode {
		"0" {
			#^ fixed constant
			if {$independentConstants == 1} {
				set meanCurrentConstant 0
				foreach previousConstantAtom $previousConstants {
					set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
					set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
				}
				set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
			} else {
				set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
				$averageFrameSelection set $constantsColumn $currentConstant
			}
		}
		"1" {
			#^ exponential increase
			if {$parametrizationCyclesCompleted > 0} {
				# This is not the first cycle
				if {$independentConstants == 1} {
					set currentConstants {}
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom*$multiplierConstant}]
						set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
						lappend currentConstants $currentConstant
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
					# Write calculated constants vector into constants column
					$averageFrameSelection set $constantsColumn $currentConstants
				} else {
					set currentConstant [expr {$previousConstant*$multiplierConstant}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			} else {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}
			
		}
		"2" {
			#^ "constant force" mode (constants adjusted such that shift by 1 A gives initially the same amount of energy
			if {$independentConstants == 1} {
				# Calculate average current constant and necessary correction
				set meanCurrentConstant 0
				foreach previousConstantAtom $previousConstants {
					set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
					set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
				}
				set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				if {$currentRmsd>0} {
					set currentConstantCorrection [expr {$energySlopeConstant/(2.0*$currentRmsd*$currentConstant)}]
				} else {
					set currentConstantCorrection 9999.999
				}
				
				# Apply correction
				set currentConstants {}
				set meanCurrentConstant 0
				foreach previousConstantAtom $previousConstants {
					set currentConstant [expr {$previousConstantAtom*$currentConstantCorrection}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					lappend currentConstants $currentConstant
					set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
				}
				set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
				# Write calculated constants vector into constants column
				$averageFrameSelection set $constantsColumn $currentConstants
			} else {
				if {$currentRmsd>0} {
					set currentConstant [expr {$energySlopeConstant/(2.0*$currentRmsd)}]
				} else {
					set currentConstant 9999.999
				}
				set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
				$averageFrameSelection set $constantsColumn $currentConstant
			}
		}
		"3" {
			#^ RMSD gradient
			set gradientMultiplierConstant [expr {($currentRmsd/$previousRmsd)>$RmsdGradient?$multiplierConstant:(1.0/$multiplierConstant)}]
			if {$parametrizationCyclesCompleted > 0} {
				# This is not the first cycle
				if {$independentConstants == 1} {
					set currentConstants {}
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom*$gradientMultiplierConstant}]
						set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
						lappend currentConstants $currentConstant
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
					# Write calculated constants vector into constants column
					$averageFrameSelection set $constantsColumn $currentConstants
				} else {
					set currentConstant [expr {$previousConstant*$gradientMultiplierConstant}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			} else {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}
			
		}
		"4" {
			#^ "oscillating constant" : spring constant amount will periodically increase and decrease exponentially (linear in log scale) between the $oscillatingConstantMin and $oscillatingConstantMax with the rate of $multiplierConstant/step and $oscillatingConstantPeriod (once in a oscillating period it will go up and down)
			# Calculate the number of steps of exponential growth/decay of the oscillating spring constant. 
			# Every $oscillatingConstantPeriod minimization cycles the ocillating constants will scale up for $oscillatingConstantSteps with $multiplierConstant constant increase rate, and then the likewise down.
			# If starting constants were different, they will remain different, but scale proportionally.
			set oscillatingConstantSteps [expr {round(log($oscillatingConstantMax/$oscillatingConstantMin)/log($multiplierConstant))}]
			set oscillatingConstant [expr {($parametrizationCyclesCompleted % $oscillatingConstantPeriod >= $oscillatingConstantSteps) ? (($parametrizationCyclesCompleted % $oscillatingConstantPeriod >= ($oscillatingConstantPeriod - $oscillatingConstantSteps)) ? $multiplierConstant : 1) : (1.0 / $multiplierConstant)}]
			
			if {$parametrizationCyclesCompleted == 0} {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}

			# Adjust the constant
			if {$independentConstants == 1} {
				set currentConstants {}
				set meanCurrentConstant 0
				foreach previousConstantAtom $previousConstants {
					set currentConstant [expr {$previousConstantAtom*$oscillatingConstant}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					lappend currentConstants $currentConstant
					set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
				}
				set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				# Write calculated constants vector into constants column
				$averageFrameSelection set $constantsColumn $currentConstants
			} else {
				set currentConstant [expr {$previousConstant*$oscillatingConstant}]
				set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
				$averageFrameSelection set $constantsColumn $currentConstant
			}
		}
		"5" {
			#^ "linear increase". Coefficient=initialConstant+multiplierConstant*timestep
			if {$parametrizationCyclesCompleted > 0} {
				# This is not the first cycle
				if {$independentConstants == 1} {
					set currentConstants {}
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
 						set currentConstant [expr {$previousConstantAtom+$multiplierConstant}]
						set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
						lappend currentConstants $currentConstant
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
					# Write calculated constants vector into constants column
					$averageFrameSelection set $constantsColumn $currentConstants
				} else {
					set currentConstant [expr {$previousConstant+$multiplierConstant}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			} else {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}
			
		}
		"6" {
			# "parabolic increase". Coefficient=initialConstant+(multiplierConstant*timestep)^2
			if {$parametrizationCyclesCompleted > 0} {
				# This is not the first cycle
				if {$independentConstants == 1} {
					set currentConstants {}
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom+pow($multiplierConstant,2)*(2*$parametrizationCyclesCompleted-1.0)}]
						set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
						lappend currentConstants $currentConstant
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
					# Write calculated constants vector into constants column
					$averageFrameSelection set $constantsColumn $currentConstants
				} else {
					set currentConstant [expr {$previousConstant+pow($multiplierConstant,2)*(2*$parametrizationCyclesCompleted-1.0)}]
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			} else {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}
			
		}
		"7" {
			# <7> - "square-hyperbolic increase" Coefficient=initialConstant/(1-multiplierConstant*timestep)^2. For "smooth" system should lead to the linear RMSD decrease mode. multiplierConstant*timestep should be < 1. multiplierConstant = (1-initialConstant/maxconstant)/timesteps
			if {$parametrizationCyclesCompleted > 0} {
				# This is not the first cycle
				if {$independentConstants == 1} {
					set currentConstants {}
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						if {[expr {double($multiplierConstant*$parametrizationCyclesCompleted)}]<1} {;	# Coefficient will not be infinity
							set currentConstant [expr {$previousConstantAtom*pow((1.0-$multiplierConstant*($parametrizationCyclesCompleted-1)),2)/pow((1.0-$multiplierConstant*$parametrizationCyclesCompleted),2)}]
						} else {
							set currentConstant 9999.999
						}
						set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
						lappend currentConstants $currentConstant
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
					# Write calculated constants vector into constants column
					$averageFrameSelection set $constantsColumn $currentConstants
				} else {
					if {[expr {double($multiplierConstant*$parametrizationCyclesCompleted)}]<1} {;	# Coefficient will not be infinity
						set currentConstant [expr {$previousConstant*pow((1.0-$multiplierConstant*($parametrizationCyclesCompleted-1)),2)/pow((1.0-$multiplierConstant*$parametrizationCyclesCompleted),2)}]
					} else {
						set currentConstant 9999.999
					}
					set currentConstant [expr {$currentConstant>9999.999?9999.999:($currentConstant<0.001?0.001:$currentConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			} else {
				# This is the very first cycle
				if {$independentConstants == 1} {
					set meanCurrentConstant 0
					foreach previousConstantAtom $previousConstants {
						set currentConstant [expr {$previousConstantAtom>9999.999?9999.999:($previousConstantAtom<0.001?0.001:$previousConstantAtom)}]
						set meanCurrentConstant [expr {$meanCurrentConstant + $currentConstant}]
					}
					set currentConstant [expr {$meanCurrentConstant/[llength $previousConstants]}]
				} else {
					set currentConstant [expr {$initialConstant>9999.999?9999.999:($initialConstant<0.001?0.001:$initialConstant)}]
					$averageFrameSelection set $constantsColumn $currentConstant
				}
			}
			
		}
	}
	
	# ___ Save spring constants
	$all writepdb $outputConstantsFileName
	# Make a copy for the log, if necessary
	if {$logOutputFiles == 1} {
		# file copy -force $outputConstantsFileName $outputConstantsFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted]
		switch $logOutputFormat {
			"pdb" {set logOutputFormatExtension "pdb"}
			"dcd" {set logOutputFormatExtension "dcd"}
			"namdbin" {set logOutputFormatExtension "coor"}
			default {set logOutputFormatExtension "pdb"}
		}
				
		animate write $logOutputFormat $outputConstantsFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted].$logOutputFormatExtension beg $averageFrame end $averageFrame
	}
	
	
	$selection frame $averageFrame
	$selection update
	# If there will be 2 selections used for target calculation?
	case $targetCalculationMode  {
		{0 1 2 3 4 9 10} {;	# There will be 1 selection (selection 2 will be ignored)
		}
		{5 6 7 8} {;	# There will be 2 selections. The same Spring constants will be calculated for joined selections 1 and 2
			$selection1 frame $averageFrame
			$selection1 update
			$selection2 frame $averageFrame
			$selection2 update
		}
	};	# If there will be 2 selections used?

	
	
	$all frame $averageFrame
	$all update
	
	# ________________ Calculate new target coordinates, if necessary
	if {$targetCalculationMode !=0} {
		# Read target coordinates from file, if necessary
		# Recall calculated coordinates for this
		$all set {x y z} $targetXYZ
		if {$preserveUnassignedConstants ==1} {
			# Read target coordinates from the existing file, calculate new coordinates, assign them and leave all old coordinates, which were not reassigned, in place
			# Useful when there are some additional harmonic restraints in the system, not related to parametryzation
			if {[catch {
				# Try to read old target coordinates
				mol new $structure type psf
				mol addfile $outputCoordinatesFileName type pdb waitfor all
			}] == 0} {
				# Old target coordinates were read
				# Remember coordinates of the selection
				set xyzCoordinates [$selection get {x y z}]
				$all set {x y z} [[atomselect top all] get {x y z}]
				$selection set {x y z} $xyzCoordinates
				puts "target coordinates read"
				mol delete top
			} else {
				# There was not constants file. Initialize one
				# Zero the whole column for the future target coordinate values. 
				if {[molinfo top get numframes]==0} {
					# unsuccessfully loaded PSF
					mol delete top
				}
				set xyzCoordinates [$selection get {x y z}]
				$all set x 0
				$all set y 0
				$all set z 0
				$selection set {x y z} $xyzCoordinates
			}
			set previousConstants [$averageFrameSelection get $constantsColumn]
		} elseif {$preserveUnassignedConstants ==0} {
			# All old sonstants are set to 0
			set xyzCoordinates [$selection get {x y z}]
			$all set x 0
			$all set y 0
			$all set z 0
			$selection set {x y z} $xyzCoordinates
		} else {
			# coordinates from the current protein snapshot will be inserted into unassigned target coordinate positions
		}
		# Calculate new target coordinates
		# Mode of the target coordinates calculation
		switch $targetCalculationMode {
			"0" {;	# <0> - "fixed target" mode: no target update
				# Do nothing (this line should never execute)
				
			}
			"1" {;	# <1> - "cylindrical projection target": target points are mapped onto cyllinder of predefined (in the "update_parameters.vmd" script) radius. They have the same angular coordinates as in the input coordinates file
				set xyCoordinates [$selection get {x y}]
				set xyCoordinatesNew {}
				
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecscale $xyAtom [expr {$cylinderRadius/[veclength $xyAtom]}]]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat [vecscale $xyAtom [expr {$cylinderRadius/[veclength $xyAtom]}]] 0]] 0 1]
					}
				}
				$selection set {x y} $xyCoordinatesNew
			}
			"2" {;	# <2> - "constant force outward": all the restrained atoms get the same amount of force in the direction from the central (z) axis through the atom.
			puts "# <2> - constant force outward"
				set xyCoordinates [$selection get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecnorm $xyAtom]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat [vecnorm $xyAtom] 0]] 0 1]
					}
				}
				$all set x 0
				$all set y 0
				$all set z 0
				$selection set {x y} $xyCoordinatesNew
			}
			"3" {;	# <3> - "Coriolis force": all the restrained atoms get the force propotional to the square of the distance to the central axis. Force will be in the direction from the central (z) axis through the atom
				# It just leaves the old x and y coordinates in place. Being miltiplied in NAMD by the spring constant (stored in the Occupancy column) it gives the Coriolis force proportional to  the radial position
				set xyCoordinates [$selection get {x y}]
				$all set x 0
				$all set y 0
				$all set z 0
				if {[catch {set twist}]} {
					# No twist
					$selection set {x y} $xyCoordinates
				} else {
					# With twist
					set xyCoordinatesNew {}
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat $xyAtom 0]] 0 1]
					}
					$selection set {x y} $xyCoordinatesNew
				}
			}
			"4" {;	# <4> - "constant shift": target coordinate of the atom shifts every time in radial direction by the fixed distance relative to the current coordinate of the atom
				set xyCoordinates [$selection get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecscale $xyAtom [expr {([veclength $xyAtom]+$constantShift)/[veclength $xyAtom]}]]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecadd $xyAtom [lrange [vectrans [transaxis z $twist] [concat [vecscale $xyAtom [expr {$constantShift/[veclength $xyAtom]}]] 0]] 0 1]]
					}
				}
				$selection set {x y} $xyCoordinatesNew
			}
			"5" {;	# <5> - "cylindrical projection (selection 1) and axis (selection 2) target": target points are mapped onto cyllinder of predefined (in the "update_parameters.vmd" script) radius. They have the same angular coordinates as in the input coordinates file
				set xyCoordinates [$selection1 get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecscale $xyAtom [expr {$cylinderRadius/[veclength $xyAtom]}]]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat [vecscale $xyAtom [expr {$cylinderRadius/[veclength $xyAtom]}]] 0]] 0 1]
					}
				}
				$selection1 set {x y} $xyCoordinatesNew
				
				if {[catch {set twist2}]} {
					# No twist
					$selection2 set x 0
					$selection2 set y 0
				} else {
					# With twist
					set xyCoordinates [$selection2 get {x y}]
					set xyCoordinatesNew {}
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecsub $xyAtom [lrange [vectrans [transaxis z [expr {-$twist2}]] [concat $xyAtom 0]] 0 1]]
					}
					$selection2 set {x y} $xyCoordinatesNew
				}
				
			}
			"6" {;	# <6> - "constant force outward (selection 1) and to the axis (selection 2)": all the restrained atoms get the same amount of force in the direction from the central (z) axis through the atom.
				set xyCoordinates [$selection1 get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecnorm $xyAtom]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat [vecnorm $xyAtom] 0]] 0 1]
					}
				}
				
				set xyCoordinates [$selection2 get {x y}]
				set xyCoordinatesNew2 {}
				if {[catch {set twist2}]} {
					# No twist - coordinates w
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew2 [vecscale [vecnorm $xyAtom] -1]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew2 [vecscale [lrange [vectrans [transaxis z [expr {-$twist2}]] [concat [vecnorm $xyAtom] 0]] 0 1] -1]
					}
				}
				
				$all set x 0
				$all set y 0
				$all set z 0
				$selection1 set {x y} $xyCoordinatesNew
				$selection2 set {x y} $xyCoordinatesNew2
			}
			"7" {;	# <7> - "Coriolis force outward (selection 1) and to the axis (selection 2)": all the restrained atoms get the force propotional to the square of the distance to the central axis. Force will be in the direction from the central (z) axis through the atom
				# It just leaves the old x and y coordinates in place. Being miltiplied in NAMD by the spring constant (stored in the Occupancy column) it gives the Coriolis force proportional to  the radial position
				set xyCoordinates [$selection1 get {x y}]
				set xyCoordinates2 [$selection2 get {x y}]
				$all set x 0
				$all set y 0
				$all set z 0
				if {[catch {set twist}]} {
					# No twist
					$selection1 set {x y} $xyCoordinates
				} else {
					# With twist
					set xyCoordinatesNew {}
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [lrange [vectrans [transaxis z $twist] [concat $xyAtom 0]] 0 1]
					}
					$selection1 set {x y} $xyCoordinatesNew
				}
				
				if {[catch {set twist2}]} {
					# No twist
					set xyCoordinatesNew2 {}
					foreach xyAtom $xyCoordinates2 {
						lappend xyCoordinatesNew2 [vecscale $xyAtom -1]
					}
				} else {
					# With twist
					set xyCoordinatesNew2 {}
					foreach xyAtom $xyCoordinates2 {
						lappend xyCoordinatesNew2 [vecscale [lrange [vectrans [transaxis z [expr {-$twist2}]] [concat $xyAtom 0]] 0 1] -1]
					}
				}
				$selection2 set {x y} $xyCoordinatesNew2
			}
			"8" {;	# <8> - "constant shift outward (selection 1) and to the axis (selection 2)": target coordinate of the atom shifts every time in radial direction by the fixed distance relative to the current coordinate of the atom
				set xyCoordinates [$selection1 get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecscale $xyAtom [expr {([veclength $xyAtom]+$constantShift)/[veclength $xyAtom]}]]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecadd $xyAtom [lrange [vectrans [transaxis z $twist] [concat [vecscale $xyAtom [expr {$constantShift/[veclength $xyAtom]}]] 0]] 0 1]]
					}
				}
				$selection1 set {x y} $xyCoordinatesNew
				
				set xyCoordinates [$selection2 get {x y}]
				set xyCoordinatesNew {}
				if {[catch {set twist2}]} {
					# No twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecscale $xyAtom [expr {([veclength $xyAtom]-$constantShift)/[veclength $xyAtom]}]]
					}
				} else {
					# With twist
					foreach xyAtom $xyCoordinates {
						lappend xyCoordinatesNew [vecsub $xyAtom [lrange [vectrans [transaxis z [expr {-$twist2}]] [concat [vecscale $xyAtom [expr {$constantShift/[veclength $xyAtom]}]] 0] ] 0 1]]
					}
				}
				$selection2 set {x y} $xyCoordinatesNew
			}
			"9" {;	# <9> - "restrained center of mass of every subselection" : atoms of every subselection provided as a sublist in the "selection" variable (for one selection use one sublist like {{resname ALA}} ) will be harmonically restrained to tranlsated location such that COM of every subselection will be restrained to it's location in the reference file

				# Read the reference structure
				if {[catch {
					mol new $structure type psf
					mol addfile $referenceCoordinates type $referenceCoordinatesFileType waitfor all
					set referenceMolID [molinfo top get id]
					puts "Reference structure loaded"
				}]} {error "Can't load the reference structure"}

				# Go trough the subselections, calculate COMs and set target coordinates
				set averageDistance 0
				set dictsncesNumber 0
				foreach subselectionText $selectionText {
					set selectionTarget [atomselect $coordinatesMolID $subselectionText frame $averageFrame]
					set selectionReference [atomselect $referenceMolID $subselectionText]
					set ComTarget [measure center $selectionTarget weight mass]
					set ComReference [measure center $selectionReference weight mass]
					set distance [vecsub $ComReference $ComTarget]
					$selectionTarget moveby $distance
					set distance [veclength $distance]
					set averageDistance [expr {$averageDistance+pow($distance,2)}]
					incr dictsncesNumber
					$selectionTarget delete
					$selectionReference delete
				}
				set currentRmsd [expr {pow($averageDistance/$dictsncesNumber,0.5)}]
# 				set selectionReference [atomselect $referenceMolID [join $selectionText " or "]]
# 				catch {set currentRmsd [measure rmsd $selection $selectionReference]}
# 				$selectionReference delete
				mol delete $referenceMolID
			}
			"10" {;	# <10> - "restrained center of mass of every fragment in selection" : atoms of every fragment in the "selection" will be harmonically restrained to tranlsated location such that COM of every fragment will be restrained to it's location in the reference file

				# Read the reference structure
				if {[catch {
					mol new $structure type psf
					mol addfile $referenceCoordinates type $referenceCoordinatesFileType waitfor all
					set referenceMolID [molinfo top get id]
					puts "Reference structure loaded"
				}]} {error "Can't load the reference structure"}
				
				# Go trough the subselections, calculate COMs and set target coordinates
				set averageDistance 0
				set dictsncesNumber 0
				set fragmentsList [lsort -unique [$selection get fragment]]
				foreach fragment $fragmentsList {
					set selectionTarget [atomselect $coordinatesMolID "($selectionText) and fragment $fragment" frame $averageFrame]
					set selectionReference [atomselect $referenceMolID "($selectionText) and fragment $fragment"]
					set ComTarget [measure center $selectionTarget weight mass]
					set ComReference [measure center $selectionReference weight mass]
					set distance [vecsub $ComReference $ComTarget]
					$selectionTarget moveby $distance
					set distance [veclength $distance]
					set averageDistance [expr {$averageDistance+pow($distance,2)}]
					incr dictsncesNumber
					$selectionTarget delete
					$selectionReference delete
				}
				set currentRmsd [expr {pow($averageDistance/$dictsncesNumber,0.5)}]
# 				set selectionReference [atomselect $referenceMolID $selectionText]
# 				catch {set currentRmsd [measure rmsd $selection $selectionReference]}
# 				$selectionReference delete
				mol delete $referenceMolID
			}
			default {
				# Unrecognized option
				puts "Unrecognized mode of target calculation"
				break
			}
		}
		# Calculated target coordinates will be savet later, after constants will be calculated (because constant force parameters require them both in  one file
		# Remember calculated coordinates for this
	}
	
	
	
		# ________ Save new target coordinates
	puts "Saving structure..."
	$all frame $averageFrame
	$all update
	# Finalize the preparation of the structure file
	case $targetCalculationMode {
		{2 3 6 7} {;	# Constant force type restraints
			# Insert spring coefficient into occupancy column. The coefficient should be in the range 0.01-999.99
			$all set occupancy 0
			$averageFrameSelection set occupancy $currentConstant
		}
	}
	
	
	if {($targetCalculationMode>0) && [catch {
			$all writepdb $outputCoordinatesFileName
			# Make a copy for the log, if necessary
			if {$logOutputFiles == 1} {
				# file copy -force $outputCoordinatesFileName $outputCoordinatesFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted]
				switch $logOutputFormat {
					"pdb" {set logOutputFormatExtension "pdb"}
					"dcd" {set logOutputFormatExtension "dcd"}
					"namdbin" {set logOutputFormatExtension "coor"}
					default {set logOutputFormatExtension "pdb"}
				}
						
				animate write $logOutputFormat $outputCoordinatesFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted].$logOutputFormatExtension beg $averageFrame end $averageFrame
			}
		}] != 0} {
		$all writepdb [molinfo top get name]tgt.pdb
	}
		
		# ___ Compose values for parametrization state file
	# Number of minimization cycles completed up to the moment
	set parametrizationCyclesCompleted $parametrizationCyclesCompleted
	# Average RMSD between the monomers and parametric targed, obtained in the last step
	set previousRmsd $currentRmsd
	# Spring constant value attained in the previous step
	set previousConstant $currentConstant
	
	# ___ Save parametrization state file
	set fileLook [open $statusFileName w]
	puts $fileLook "# Autogenerated minimization status file"
	puts $fileLook "parametrizationCyclesCompleted\t$parametrizationCyclesCompleted"
	puts $fileLook "previousRmsd\t$previousRmsd"
	puts $fileLook "previousConstant\t$previousConstant"
	close $fileLook
	# Make a copy for the log, if necessary
	if {$logOutputFiles == 1} {
# 		file copy -force $statusFileName $statusFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted]
		set fileLook [open $statusLogFileName a]
		puts $fileLook "$parametrizationCyclesCompleted\t$previousRmsd\t$previousConstant"
		close $fileLook
	} 
	
				
	# _____ Check for the exit imperative
	if {$parametrizationCyclesCompleted >= $parametrizationStepsMax} {
		# Maximal number of parametrization steps reached
		puts "Maximal number of parametrization steps reached"
		file copy -force $outputCoordinatesFileName [join "$outputCoordinatesFileName _final.pdb" {}]
		
	} else {
		if {$currentRmsd < $targetRmsd} {
			# Target RMSD from the parametric structure reached
			switch $constantsScaleMode {
				"0" -
				"1" -
				"2" -
				"3" {
					puts "Current RMSD $currentRmsd:   Target RMSD ($targetRmsd) reached"
					file copy -force $outputCoordinatesFileName [join "$outputCoordinatesFileName _final.pdb" {}]
				}
			}
		}
	}
} else {
	# No parametrization was performed
	
	# ___ Save parametrization state file
	set fileLook [open $statusFileName w]
	puts $fileLook "# Autogenerated minimization status file"
	puts $fileLook "parametrizationCyclesCompleted\t$parametrizationCyclesCompleted"
	puts $fileLook "previousRmsd\t0"
	puts $fileLook "previousConstant\t0"
	close $fileLook
	# Make a copy for the log, if necessary
	if {$logOutputFiles == 1} {
# 		file copy -force $statusFileName $statusFileName[format "_%0[string length $parametrizationStepsMax]i" $parametrizationCyclesCompleted]
		set fileLook [open $statusLogFileName a]
		puts $fileLook "$parametrizationCyclesCompleted\t0\t0"
		close $fileLook
	} 
				
	# _____ Check for the exit imperative
	if {$parametrizationCyclesCompleted >= $parametrizationStepsMax} {
		# Maximal number of parametrization steps reached
		puts "Maximal number of parametrization steps reached"
		file copy -force $outputCoordinatesFileName [join "$outputCoordinatesFileName _final.pdb" {}]
	}
}

cd $oldFolder
#mol delete all

puts "Finished!!!"

set spent_time [expr {[clock seconds]-$start_time}];
set spent_time_h [expr {int(floor($spent_time/3600))}];
set spent_time_m [expr {int(floor(($spent_time-$spent_time_h*3600)/60))}];
set spent_time_s [expr {int($spent_time-$spent_time_h*3600-$spent_time_m*60)}];

puts "Time spent:  $spent_time_h hours  $spent_time_m min  $spent_time_s sec"

exit