rattle.tcl

MD中rattle 算法

proc rattle {type args} {
#Initial check of type value before other errors can occur
    if {$type != "off" && $type != "bonds" && $type != "angles" && $type != "water"} {
	error "rattle: unknown type \"$type\": must be `bonds', `angles', `water' or `off'"
    }


    #set system_h
    database handle system_h System.

    #reset the resetEnergyExpression and chargeChanged flags under the
    #Global table in the system database
    System_SetGlobal resetEnergyExpression 1
    System_SetGlobal chargeChanged 1

    set rattleMsg ""
    vector bl {}
    vector av {}
#Check/perform the `off' function
    database handle dbh {} dbn
    object dbn cast int

    $l = [llength $args]
    if {$type == "off"} {
	if $l {
	    error "rattle off: extra arguments given"
	}
	rattleSystem off [object dbn]
	return
    }
    rattleSystem on [object dbn]

#Check/grab the tolerance arguments - the optional values are considered further on
    if {[lindex $args 0] == "-tolerance"} {
	if {$l < 2} {
	    error "rattle: with -tolerance: no tolerance value specified"
	}
	vector tol [lindex $args 1]
	if {[object tol info type] != "double"} {
	    error "rattle: tolerance value must be a vector of reals"
	}
	incr l -2
	set args [lrange $args 2 end]
    } else {
	vector tol 1.0e-7
    }
#Check/set up water parameters
    if {$type == "water"} {
    #Check/set up water ff type
	global env
	if {[lindex $args 0] == "-waterType"} {
	    if {$l < 2} {
		error "rattle water: option value not specified for `-waterType'"
	    }
	    set wtype [string toupper [lindex $args 1]]
	    incr l -2
	    set args [lrange $args 2 end]
	} else {
	    set wtype SPC
	}
	if {[regexp -nocase esff $env(FORCEFIELD)]} {
	    if {$wtype != "CURRENT"} {
		error "rattle: ESFF forcefield does not support SPC or TIP3P waters"
	    }
	}
    #Check for tolerance again in case this comes now
	if {[lindex $args 0] == "-tolerance"} {
	    if {$l < 2} {
		error "rattle: with -tolerance: no tolerance value specified"
	    }
	    vector tol [lindex $args 1]
	    if {[object tol info type] != "double"} {
		error "rattle: tolerance value must be a vector of reals"
	    }
	    incr l -2
	    set args [lrange $args 2 end]
	}
    #Check no more arguments provided - waters is all or nothing
	if $l {
	    error "rattle water: too many arguments provided: no water subsets may be specified"
	}
    #Grab target atom sets - must be waters
	$dbh select -nosubdivide {WTR OHH} Monomer.Type waters
	if ![object waters info nItems] {
	    echo "Warning: rattle: system does not contain water molecules (monomer type \"WTR\" or \"OHH\")"
	    return
	}
	$dbh select -nosubdivide $waters Atom.Monomer waters
	$dbh select -nosubdivide -regexp ^O.*   Atom.Name waterO  $waters
	$dbh select -nosubdivide -regexp ^H.*1$ Atom.Name waterH1 $waters
	$dbh select -nosubdivide -regexp ^H.*2$ Atom.Name waterH2 $waters
    #Check/setup water ff parameters and geometries
	if {$wtype == "SPC"} {
	    vector bl 1.00
	    vector wc -0.820
	    molGeom set distance $bl $waterO $waterH1
	    molGeom set distance $bl $waterO $waterH2
	    molGeom set angle 109.47 $waterH1 $waterO $waterH2
	    $dbh set ospc Atom.Type $waterO
	    $dbh set hspc Atom.Type $waterH1
	    $dbh set hspc Atom.Type $waterH2
	    $dbh set $wc Atom.Charge $waterO
	    $dbh set 0.410 Atom.Charge $waterH1
	    $dbh set 0.410 Atom.Charge $waterH2
	} elseif {$wtype == "TIP3P"} {
	    vector bl 0.957
	    vector wc -0.834
	    molGeom set distance $bl $waterO $waterH1
	    molGeom set distance $bl $waterO $waterH2
	    molGeom set angle 104.52 $waterH1 $waterO $waterH2
	    $dbh set otip Atom.Type $waterO
	    $dbh set htip Atom.Type $waterH1
	    $dbh set htip Atom.Type $waterH2
	    $dbh set $wc Atom.Charge $waterO
	    $dbh set 0.417 Atom.Charge $waterH1
	    $dbh set 0.417 Atom.Charge $waterH2
	} elseif {$wtype != "CURRENT"} {
	    error "rattle water: unknown forcefield specification \"$wtype\": must be SPC, TIP3P or CURRENT"
	}
    #Put OH distance into the Rattle table
	select atoms1 $waterH1
	object atoms1 append $waterH2
	select atoms2 $waterO
	object atoms2 append $waterO
	if {$wtype == "CURRENT"} {
	    molGeom get distance bl $atoms1 $atoms2
	}
	addRattle $atoms1 $atoms2 $bl $tol
    #Set up Atom lists for putting into the Rattle table for atoms adjacent to angle
	set atoms1 $waterH1
	set atoms2 $waterH2
	vector bl ""
    } elseif {$type == "bonds" } {
#Check/set up bond values and atom sets
    #Check for value specification
	if {[lindex $args 0] == "-bondLength"} {
	    if {$l < 2} {
		error "rattle bonds: option value not specified for `-bondLength'"
	    }
	    vector bl [lindex $args 1]
	    incr l -2
	    set args [lrange $args 2 end]
	}
    #Check for tolerance again in case this comes now
	if {[lindex $args 0] == "-tolerance"} {
	    if {$l < 2} {
		error "rattle: with -tolerance: no tolerance value specified"
	    }
	    vector tol [lindex $args 1]
	    if {[object tol info type] != "double"} {
		error "rattle: tolerance value must be a vector of reals"
	    }
	    incr l -2
	    set args [lrange $args 2 end]
	}
	switch $l {
	0 {
	    $dbh select Water Molecule.Type waters
	    if [object waters info nItems] {
		select nWater "!(*<Type>Water):*:Atom;*"
		if [object nWater info nItems] {
		    getBondAtoms atoms1 atoms2 $nWater
		    set rattleMsg "constrained all bonds for all non-water molecules"
		} else {
		    getBondAtoms atoms1 atoms2
		    set rattleMsg "constrained all bonds for all water molecules"
		}
	    } else {
		getBondAtoms atoms1 atoms2
		set rattleMsg "constrained all bonds for all molecules"
	    }
	  }
	1 {
	    set subtype [string tolower [subset get subs [lindex $args 0]]]
	    if {$subtype != "distance"} {
		error "rattle bonds: with 1 subset specification: must be a `distance' subset"
	    }
	    $dbh get atoms1 Distance.Point1 $subs
	    $dbh get atoms2 Distance.Point2 $subs
	  }
	2 {
	    select atoms1 [lindex $args 0] ; $natoms1 = [object atoms1 info nItems]
	    select atoms2 [lindex $args 1] ; $natoms2 = [object atoms2 info nItems]
	    if {$natoms1 != $natoms2} {
		error "rattle bonds: atom lists must specify the same number of items"
	    }
	  }
	default {error "rattle bonds: too many atom specifications provided"}
	}
    } elseif {$type == "angles"} {
#Check/set up angle values and atom sets
    #Check for value specification
	if {[lindex $args 0] == "-angleValue"} {
	    if {$l < 2} {
		error "rattle angles: option value not specified for `-angleValue'"
	    }
	    vector av [lindex $args 1]
	    incr l -2
	    set args [lrange $args 2 end]
	}
    #Check for tolerance again in case this comes now
	if {[lindex $args 0] == "-tolerance"} {
	    if {$l < 2} {
		error "rattle: with -tolerance: no tolerance value specified"
	    }
	    vector tol [lindex $args 1]
	    if {[object tol info type] != "double"} {
		error "rattle: tolerance value must be a vector of reals"
	    }
	    incr l -2
	    set args [lrange $args 2 end]
	}
	switch $l {
	0 {
	    $dbh select Water Molecule.Type waters
	    if [object waters info nItems] {
		select nWater "!(*<Type>Water):*:Atom;*"
		if [object nWater info nItems] {
		    getAngleAtoms atoms1 atoms2 atoms3 $nWater
		    set rattleMsg "constrained all angles for all non-water molecules"
		} else {
		    getAngleAtoms atoms1 atoms2 atoms3
		    set rattleMsg "constrained all angles for all water molecules"
		}
	    } else {
		getAngleAtoms atoms1 atoms2 atoms3
		set rattleMsg "constrained all angles for all molecules"
	    }
	  }
	1 {
	    set subtype [string tolower [subset get subs [lindex $args 0]]]
	    if {$subtype != "angle"} {
		error "rattle angles: with 1 subset specification: must be an `angle' subset"
	    }
	    $dbh get atoms1 Angle.Point1 $subs
	    $dbh get atoms2 Angle.Point2 $subs
	    $dbh get atoms3 Angle.Point3 $subs
	  }
	3 {
	    select atoms1 [lindex $args 0] ; $natoms1 = [object atoms1 info nItems]
	    select atoms2 [lindex $args 1] ; $natoms2 = [object atoms2 info nItems]
	    select atoms3 [lindex $args 2] ; $natoms3 = [object atoms3 info nItems]
	    if {$natoms1 != $natoms2 || $natoms1 != $natoms3} {
		error "rattle angles: atom lists must specify the same number of items"
	    }
	  }
	default {error "rattle angles: incorrect number of atom specifications provided"}
	}
    #Fix specified angle values
	if [object av info nItems] {
	    if [catch {molGeom set angle $av $atoms1 $atoms2 $atoms3} msg] {
		set msg [string range $msg 14 end]
		error "rattle angles: $msg"
	    }
	}
	set atoms2 $atoms3
    }
#Fix/get specified bonds lengths
    if {[object bl info nItems] && $type == "bonds"} {
	if [catch {molGeom set distance $bl $atoms1 $atoms2} msg] {
	    if {[string first "directly bonded" $msg] || [string first "ring structure" $msg]} {
		set msg [string range $msg 14 end]
		echo Warning: rattle bonds: unable to set geometry: $msg
	    } else {
		error "rattle bonds: $msg"
	    }
	}
    } else {
	molGeom get distance bl $atoms1 $atoms2
    }
#Finally fill in the Rattle table
    addRattle $atoms1 $atoms2 $bl $tol
    if {$rattleMsg != ""} {echo "rattle: $rattleMsg"}
}
proc addRattle {atms1 atms2 bondLength tolerance} {
#This procedure assumes that atoms1 & atoms2 are atom objects
#and that bondLength and tolerance are scalar vectors
#Select (copy) the atom lists and check/create the Rattle table
    select atoms1 $atms1   ; $nat1 = [object atoms1 info nItems]
    select atoms2 $atms2
    vector bdl $bondLength ; $nbdl = [object bdl info nItems]
    vector tol $tolerance  ; $ntol = [object tol info nItems]
    if {$nbdl != 1 && $nbdl != $nat1} {
	error "rattle: constraint value must contain 1 or $nat1 values"
    }
    if {$ntol != 1 && $ntol != $nat1} {
	error "rattle: torelance value must contain 1 or $nat1 values"
    }
    database handle dbh $atoms1
    if !{[$dbh exists Rattle]} {
	$dbh createTable Rattle
	$dbh createColumn -nodeleteref Rattle.Atom-1 reference Atom
	$dbh createColumn -nodeleteref Rattle.Atom-2 reference Atom
	$dbh createColumn Rattle.BondLength double 
	$dbh createColumn Rattle.Tolerance double 1.0e-7
	$dbh create [object atoms1 info nItems] Rattle newRows
	$dbh set $atoms1 Rattle.Atom-1 $newRows
	$dbh set $atoms2 Rattle.Atom-2 $newRows
	$dbh set $bdl Rattle.BondLength $newRows
	$dbh set $tol Rattle.Tolerance $newRows
	return
    }
#Filter out all the atom pairs already in the table...
#First get all rows for which atoms1 & atoms2 pairs are in
    $dbh select $atoms1 Rattle.Atom-1 a1
    $dbh select -itempergroup $atoms2 Rattle.Atom-2 a2 $a1
    if [object a2 info nItems] {
    #Use the grouping on a2 to filter the initial atom-pair lists
	object a2 group -size f4
    #Apply filters to bondLength and fill in for old row values
	if {[object bdl info nItems] == 1} {
	    $dbh set $bdl Rattle.BondLength $a2
	} else {
	    object bl get $bdl
	    object bl filter $f4
	    $dbh set $bl Rattle.BondLength $a2
	}
    #Apply filters to tolerance and fill in for old row values
	if {[object tol info nItems] == 1} {
	    $dbh set $tol Rattle.Tolerance $a2
	} else {
	    object tv get $tol
	    object tv filter $f4
	    $dbh set $tv Rattle.Tolerance $a2
	}
    #Get all new atom pairs
	object f4 cast boolean
	vector f5 not $f4
	object atoms1 filter $f5
	object atoms2 filter $f5
    #Adjust tol and bdl
	if {[object bdl info nItems] > 1} {
	    object bdl filter $f5
	}
	if {[object tol info nItems] > 1} {
	    object tol filter $f5
	}
    }
#Filter out all the atom pairs already in the table in the oposite order
    if [object atoms2 info nItems] {
	$dbh select $atoms2 Rattle.Atom-1 a1
	$dbh select -itempergroup $atoms1 Rattle.Atom-2 a2 $a1
    } else {
	vector a2 ""
    }
    if [object a2 info nItems] {
    #Use the grouping on a2 to filter the initial atom-pair lists
	object a2 group -size f4
    #Apply filters to bondLength and fill in for old row values
	if {[object bdl info nItems] == 1} {
	    $dbh set $bdl Rattle.BondLength $a2
	} else {
	    object bl get $bdl
	    object bl filter $f4
	    $dbh set $bl Rattle.BondLength $a2
	}
    #Apply filters to tolerance and fill in for old row values
	if {[object tol info nItems] == 1} {
	    $dbh set $tol Rattle.Tolerance $a2
	} else {
	    object tv get $tol
	    object tv filter $f4
	    $dbh set $tv Rattle.Tolerance $a2
	}
    #Get all new atom pairs
	object f4 cast boolean
	vector f5 not $f4
	object atoms1 filter $f5
	object atoms2 filter $f5
    #Adjust tol and bdl
	if {[object bdl info nItems] > 1} {
	    object bdl filter $f5
	}
	if {[object tol info nItems] > 1} {
	    object tol filter $f5
	}
    }
#Create the new Rattle rows and fill in
    if [object atoms1 info nItems] {
	$dbh create [object atoms1 info nItems] Rattle newRows
	$dbh set $atoms1 Rattle.Atom-1 $newRows
	$dbh set $atoms2 Rattle.Atom-2 $newRows
	$dbh set $bdl Rattle.BondLength $newRows
	$dbh set $tol Rattle.Tolerance $newRows
    }
}
proc rattleSystem {op dbn} {
    global RattleSystemData
    if [catch {set RattleSystemData($dbn)}] {
	if {$op == "off"} {return} elseif {$op != "on"} {
	    error "rattleSystem: arg#1 <op>: operation must be `on' or `off'"
	}
	database handle dbh #$dbn
    #Mark Molecule.Type column for all water molecules
	$dbh select -nosubdivide {WTR OHH} Monomer.Type waters
	if [object waters info nItems] {
	    $dbh get water Monomer.Molecule $waters
	    $dbh set Water Molecule.Type $water
	    object water range $waters 0
	    $dbh select $water Atom.Monomer watms
	} else {
	    set RattleSystemData($dbn) 1
	    return
	}
    #Store to global array of database original Atom charge and type values for water atoms
	$dbh get ch Atom.Charge $watms
	$dbh get at Atom.Type $watms
	set RattleSystemData($dbn) "[object ch] [object at]"
    } elseif {$op == "on"} {return} elseif {$op == "off"} {
    #Repair atom table to original form
	set x $RattleSystemData($dbn)
	database handle dbh #$dbn
	if ![catch {select aO "*<Type>Water:*:O*"}] {
	    select aH1 "*<Type>Water:*:H*1"
	    select aH2 "*<Type>Water:*:H*2"
	    $dbh set [lindex $x 0] Atom.Charge $aO
	    $dbh set [lindex $x 1] Atom.Charge $aH1
	    $dbh set [lindex $x 2] Atom.Charge $aH2
	    $dbh set [lindex $x 3] Atom.Type $aO
	    $dbh set [lindex $x 4] Atom.Type $aH1
	    $dbh set [lindex $x 5] Atom.Type $aH2
	}
    #Unset the database encountered variable (leaving Molecule as is)
	unset RattleSystemData($dbn)
	$dbh delete Rattle
    } else {
	error "rattleSystem: arg#1 <op>: operation must be `on' or `off'"
    }
}