sybdb

来自「最经典的分子对结软件」· 代码 · 共 458 行

#!/bin/sh

# ///////////////////////////////////////////////////////////////////////////
# Script to take cleaned up multi-MOL2 format file from sdf2mol2 
# and add hydrogens, compute charges, and more.  DAG 2/95
#
# Features:
#	o removal of all but the largest covalently bonded structure
#	o hydrogen addition
#		- also corrects defects in Sybyl hydrogen addition
#			e.g. R3P becomes R3PH - corrected by this script
#	o charge computation
#		- uses Gasteiger Marsili charges, including PI charges
#		- formal charge adjustments:
#			DIRECT:
#			- charge equalization of carboxylate oxygens
#			- charge equalization of sulf(on)ate oxygens
#			- charge equalization of phos(on)ate oxygens
#			- net neutral isocyanates
#			INDIRECT:
#			- charge equalization of amidinium nitrogens
#			- charge equalization of guanidinium nitrogens
#			- charge equalization of nitro oxygens
#	o renames all atoms sequentially
#
# This script operates by creating a controlling sybyl command file which
# runs an SPL macro (also generated by this script) to perform conversion.
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

# Please customize the following appropriate to your site
#----------------------------------------------------------------------------
sybcommand="/bert/trigo/trigo sybyl6.3"
TA_ROOT=/bert/sybyl63
TA_DEMO=$TA_ROOT/demo
#----------------------------------------------------------------------------


if [ $# -ne 2 ]; then
	echo "Usage:  `basename $0` <input_filename> <output_filename>"
	echo "   where <input_filename> is a MOL2 format file resulting from"
	echo "             the program sdf2mol2, and"
	echo "         <output_filename> is the desired output name for the"
	echo "             MOL2 format with hydrogens and charges."
	echo "             This output would be suitable for input to mol2db"
	echo "             for generation of a dock database."
	exit 1
fi

# set up file names
sybscript=sybdb.script
sybmacro=sybdb.spl
syboutput=sybdb.out
logfile=sybdb.log
vers=1.00

sybprefix="##sybdb##"
mol2file=$1
mol2out=$2

# Check setup
#----------------------------------------------------------------------------

if [ ! -f $mol2file ]; then
	echo "Can not find input file: $mol2file"
	exit 1
fi

if [ -f $mol2out ]; then /bin/rm -f $mol2out; fi

if [ ! -f gastpar.tab ]; then
	echo "WARNING - no gastpar.tab Gasteiger Marsili parameter file"
	echo "was found in this directory.  Unless you have modified the"
	echo "Tripos-distributed version of this file in"
	echo "       \$TA_ROOT/sybylbase/tables/gastpar.tab"
	echo "sulfoxides and sulfones will be assigned zero charges!"
	echo
fi

# create the sybyl command script
#----------------------------------------------------------------------------
echo Creating Sybyl command script...

cat > $sybscript << EndOfScript

# Add some bogus bondings so these hydrogens are not really added.  Sybyl
# screws up by adding hydrogens to N.1 and tertiary N.ar.
parameter add bond_type N.ar H 1 no | |
echo Added N.ar-H bond_type
parameter add bond_length N.ar H 1 1.0 | |
echo Added N.ar-H bond_length

parameter add bond_type N.1 H 1 no | |
echo Added N.1-H bond_type
parameter add bond_length N.1 H 1 1.0 | |
echo Added N.1-H bond_length

parameter add bond_type O.2 H 1 no | |
echo Added O.2-H bond_type
parameter add bond_length O.2 H 1 1.0 | |
echo Added O.2-H bond_length

# Set up to use Gasteiger-Marsili pi charges, which are off by default.
tailor set charge pi_iterations 6 pi_damping 0.5 | |
echo Pi charges turned on

# Load metal parameters so that metals are recognized upon read-in.
parameter open $TA_DEMO/metals.tpd |
echo Metal parameters loaded

# Load the conversion SPL macro into sybyls command language.
uims load $sybmacro

# Execute the conversion
$sybmacro

quit
EndOfScript


# Create the sybyl macro - note the use of backslashes to avoid substitution
# of shell variables for sybyl script variables.  However, do NOT escape the
# sybprefix variable so it will get substituted by sh.
#----------------------------------------------------------------------------
echo Creating Sybyl SPL macro...
cat > $sybmacro << EndOfMacro

@MACRO

$sybmacro sybylbasic y

setvar prefix $sybprefix

# read the structures in 
mol2 mult_in m1 $mol2file

#uims verify on

for area in %mols(*)

#	______________________________________________________________________
#	Set up for processing of this molecule
#	______________________________________________________________________
	setvar molname %mol_info(\$area name)
	echo \$area \$molname
	default \$area

#	______________________________________________________________________
#	Remove all but the largest substructure
#	______________________________________________________________________
	setvar nsub %mol_info(\$area nsubsts)

	if %gt(\$nsub 1)
#		check the size of each substructure - keep only the largest
		setvar maxsize 0
		setvar biggest 0
		for sub in %substs(\$area)
			setvar natm %count(%conn_atoms(%subst_info(\$sub root_atom)))
			if %gt(\$natm \$maxsize)
				setvar maxsize \$natm
				setvar biggest \$sub
			endif
		endfor

#		catenate the string of substructures to be deleted - use
#		this method instead of sequential deletion b/c the substructure
#		numbers get modified upon every deletion
		setvar subnum 0
		setvar delsubs %cat(\$area "(")
		for sub in %substs(\$area)
			if %not(%streql(\$sub \$biggest))
				setvar subnum %math(\$subnum + 1)
				if %not(%eq(\$subnum 1))
					setvar delsubs %cat("\$delsubs" "+")
				endif
				setvar delsubs %cat("\$delsubs" \$sub)
			endif
		endfor
		setvar delsubs %cat("\$delsubs" ")")
		remove substructure \$delsubs
	endif
	

#	______________________________________________________________________
#	Add hydrogens.
#	______________________________________________________________________
	fillvalence \$area H


#	______________________________________________________________________
#	Remove dummy atoms.
#	______________________________________________________________________
	for dum in %atoms(<Du>)
		remove atom \$dum
	endfor


#	______________________________________________________________________
#	Rename atoms sequentially - to make sure added hydrogens have names.
#	This is redone again later after all processing is performed.
#	______________________________________________________________________
	modify atom name \$area sequential_auto

#	Create a few sets - all atoms, hydrogen atoms, and heavy atoms
	setvar all_set %set_create(%atoms(*))
	setvar hyd_set %set_create(%atoms(<H>))
	setvar hvy_set %set_create(%set_diff(\$all_set \$hyd_set))

#	Variable nmod contains the number of atoms whos formal charge needs to
#	be modified.
	setvar nmod 0


#	______________________________________________________________________
#	Check for functionalities for which sybyl inappropriately adds Hs.
#	Remove these hydrogens:
#	______________________________________________________________________

#		Isocyanate groups should have a -1 charge on the carbon
#		to accompany the +1 on the nitrogen

	for a in %atoms(<C.1>)

		setvar atyp %atom_info(\$a type)

#		Create a set of this atoms neighbors
		setvar nset %set_create(%atom_info(\$a neighbors))

#		see if this atom has any heavy-atom neighbors
		setvar heavy_neighbors %set_and(\$hvy_set \$nset)
		if %neq(%count(\$heavy_neighbors) 0)
			continue
		endif

#		see if this carbon has 1 non-hydrogen connection - if 
#		not, skip it.
		setvar heavy_neighbors %atoms(\$heavy_neighbors)
		if %neq(%count(\$heavy_neighbors) 1)
			continue
		endif

#		see if this one heavy atom connection is indeed an sp nitrogen
#		- if not, skip it
		if %not(%streql(%atom_info(\$heavy_neighbors type) N.1))
			continue
		endif

		setvar nremove 0

#		Loop over the intersection of neighbors of this atom
#		and the set of all hydrogens = this atoms hydrogen.
		setvar hset %set_and(\$nset \$hyd_set)
		if %eq(%count(\$hset) 0)
			continue
		endif

		for aa in %atoms(\$hset)
			echo \$prefix Stripped bad H on atom \$a \$atyp
			remove atom \$aa
			setvar nremove %math(\$nremove + 1)
		endfor

#		No need to modify the formal charge as sybyl will recognize
#		that a C.1 carbon with only one neighbor should be -1.

	endfor


#		Phosphorous with 3 non-hydrogen connections should not get
#		an additional hydrogen.

	for a in %atoms(<P.3>)

		setvar atyp %atom_info(\$a type)

#		Create a set of this atoms neighbors
		setvar nset %set_create(%atom_info(\$a neighbors))

#		see if this phosphorous has 3 non-hydrogen connections - if 
#		not, skip it.  This block is designed to remove erroneously
#		added hydrogens on the phosphorous in things like tri-alkyl
#		phosphines.
		setvar hset %set_and(\$hvy_set \$nset)
		if %eq(%count(\$hset) 0)
			continue
		endif

		if %neq(%count(%atoms(\$hset)) 3)
			continue
		endif

		setvar nremove 0

#		Loop over the intersection of neighbors of this atom
#		and the set of all hydrogens = this atoms hydrogens.
		setvar hset %set_and(\$nset \$hyd_set)
		if %eq(%count(\$hset) 0)
			continue
		endif

		for aa in %atoms(\$hset)
			echo \$prefix Stripped bad H on atom \$a \$atyp
			remove atom \$aa
			setvar nremove %math(\$nremove + 1)
		endfor

#		if some hydrogens were removed, adjust the formal
#		charge on the parent atom appropriately
		if %gt(\$nremove 0)
			setvar newchg 0.0
			setvar nmod %math(\$nmod + 1)
			setvar modatm[\$nmod] \$a
			setvar modchg[\$nmod] \$newchg
		endif
	endfor



#	______________________________________________________________________
#	Now adjust formal charges - look for particular functionalities
#	______________________________________________________________________

#	Create a set of all O.co2 atoms
	setvar oco2_set %set_create(%atoms(<O.co2>))

#	check if there were any carboxylate-type oxygens
	if %eq(%count(\$oco2_set) 0)
		setvar noco2 0
	else
		setvar noco2 %count(%atoms(\$oco2_set))
	endif

#	Sulf(on)ates and Phosph(on)ates:
#	______________________________________________________________________
	if %gt(\$noco2 0)
		for a in %atoms(<S.3>,<P.3>)

			setvar atyp %atom_info(\$a type)

#			Create a set of this atoms neighbors
			setvar nset %set_create(%atom_info(\$a neighbors))

			setvar nox 0
			setvar startat \$nmod

#			Loop over the intersection of neighbors of this atom
#			and the set of all O.co2 atoms = O.co2 on this atom
			setvar oset %set_and(\$nset \$oco2_set)
			if %eq(%count(\$oset) 0)
				continue
			endif

			for aa in %atoms(\$oset)
				setvar nox %math(\$nox + 1)
				setvar nmod %math(\$nmod + 1)
				setvar modatm[\$nmod] \$aa
			endfor

#			If any O.co2s were found, adjust the charges
#			The following is a wierd way of determining
#			what the charge on the oxygens should be:
#			Phosphates and phosphonates have 3 O.co2s
#			with a net charge of -2, giving each a -0.667
#			formal charge.  Isolated phosphates have 4
#			O.co2s with a net charge of -3 -> fc of -0.75.
#			Correspondingly, sulfates and sulfonates have
#			3 O.co2s with a net charge of -1, giving each
#			a fc of -0.333, while isolated sulfate oxygens
#			should each have a charge of -0.5 (-2/4).
#			So a rule emerges (limited, perhaps):
#			   oxygen charge = (maxchg-#oxygens)/#oxygens
#				where maxchg = 1 for phosphorous
#					   and 2 for sulfur
			switch \$atyp
				case P.3)
					setvar maxchg 1
					;;
				case S.3)
					setvar maxchg 2
					;;
				case )
					setvar maxchg 0
					;;
			endswitch
			setvar oxcharge %math( (\$maxchg - \$nox) / \$nox)

#			adjust the charges of the oxygens
			setvar lo %math(\$startat + 1)
			setvar hi %math(\$startat + \$nox)
			for o in %range(\$lo \$hi)
				setvar modchg[\$o] \$oxcharge
			endfor
	
#			set the charge of the sulfur/phosphorous to 0
			setvar nmod %math(\$nmod + 1)
			setvar modatm[\$nmod] \$a
			setvar modchg[\$nmod] 0.0
		endfor
	endif

#	Catenate the string of atoms whos formal charge must be modified.
#	______________________________________________________________________
	setvar modchgstr ""
	if %gt(\$nmod 0)
		for i in %range(1 \$nmod)
			setvar a \$modatm[\$i]
			setvar modchgstr %eval("\$modchgstr \$a \$modchg[\$i]")
		endfor
	endif


#	______________________________________________________________________
#	Compute charges with modified formal charges, if any.
#	______________________________________________________________________
	charge \$area compute gasteiger \$modchgstr |

#	echo back to the user the changes in formal charges that were made
	if %gt(\$nmod 0)
		for i in %range(1 \$nmod)
			setvar a \$modatm[\$i]
			echo \$prefix Modified Formal Charge: \$i \$a %atom_info(\$a name) %atom_info(\$a type) \$modchg[\$i]
		endfor
	endif

#	______________________________________________________________________
#	Rename atoms sequentially.
#	______________________________________________________________________
	modify atom name \$area sequential_auto


endfor

#uims verify off

# write the output
mol2 mult_out m* $mol2out
.

EndOfMacro



# Run the conversion inside sybyl
#----------------------------------------------------------------------------
echo Executing conversion...

# set up to use no graphics
set TA_AUTO_TERMINAL=NOGRAPHICS

$sybcommand < $sybscript > $syboutput

echo Compiling results...
egrep '^M[0123456789]* |^   <|^Warning|^Unknown|^'$sybprefix $syboutput > $logfile

echo Done.  Please check the logfile $logfile for conversion report.


.