#!/bin/sh
#
#  This script reads in a pdb file, and generates constraints
#    for the nmr module that force the proper chirality on 
#    tetrahedral carbons.  The pdb file must have the proper
#    AMBER atom order, e.g. have been written by LEaP or ambpdb.
#
#  It also generates "trans" constraints for all peptide bonds.
#  NOTE!!!  You must edit the output of this program manually to
#     change "trans" peptide constraints to "cis" peptide constraints
#     if needed!
#
#  Usage:  makeCHIR_RST pdb-file output-constraint-file
#
if [ $# = 0 ] ; then
    echo "usage:  $0  pdb_file  output_constraint_file "
    exit 1
fi
#
#  choose an awk cmd - some machines still have the 'old'
#   awk (10 years outdated), so to be conservative,
#   see if the more recent nawk or gawk versions exist
#   and use them.
#
AWKS="nawk gawk awk"
AWK="none"
P=`echo $PATH | sed 's/:/ /g'`
for j in $AWKS ; do
    for i in $P ; do
        if [ -x ${i}/${j} ]  ; then
            AWK=$i/$j
            break
        fi
    done
    if [ $AWK != "none" ] ; then
        break
    fi
done
if [ $AWK = "none" ] ; then
    echo "$0: no awk version found in path (tried: $AWKS)"
    echo PATH = $PATH
    exit 1
fi
echo "  (using $AWK)"

cat <<eof > chiral_defs
#
# chiral definitions for PDB proton nonemclature
# the triple product for these tetrads should always be positive
#
ALA N C HA CB
#
ARG CA CG HB2 HB3
ARG CB CD HG2 HG3
ARG CG NE HD2 HD3
ARG N C HA CB
#
ASN CA CG HB2 HB3
ASN N C HA CB
#
ASP CA CG HB2 HB3
ASP N C HA CB
#
CYS CA SG HB2 HB3
CYS N C HA CB
#
CYX CA SG HB2 HB3
CYX N C HA CB
#
GLN CA CG HB2 HB3
GLN CB CD HG2 HG3
GLN N C HA CB
#
GLU CA CG HB2 HB3
GLU CB CD HG2 HG3
GLU N C HA CB
#
GLY N C HA2 HA3
#
HIP  CA CG HB2 HB3
HIP  N C HA CB
#
HID  CA CG HB2 HB3
HID  N C HA CB
#
HIE  CA CG HB2 HB3
HIE  N C HA CB
#
HIS  CA CG HB2 HB3
HIS  N C HA CB
#
ILE CA CG2 CG1 HB
ILE CB CD1 HG12 HG13
ILE N C HA CB
#
LEU CA CG HB2 HB3
LEU CB CD1 CD2 HG
LEU N C HA CB
#
LYS CA CG HB2 HB3
LYS CB CD HG2 HG3
LYS CD NZ HE2 HE3
LYS CG CE HD2 HD3
LYS N C HA CB
#
MET CA CG HB2 HB3
MET CB SD HG2 HG3
MET N C HA CB
#
PHE CA CG HB2 HB3
PHE N C HA CB
#
PRO CA CG HB2 HB3
PRO CB CD HG2 HG3
PRO CG N  HD2 HD3
PRO N C HA CB
#
SER CA OG HB2 HB3
SER N C HA CB
#
THR CA CG2 OG1 HB
THR N C HA CB
#
TRP CA CG HB2 HB3
TRP N C HA CB
#
TYR CA CG HB2 HB3
TYR N C HA CB
#
VAL CA CG1 CG2 HB
VAL N C HA CB
#
DC  C2' O4' N1 H1'
DC  C3' C1' H2'1 H2'2
DC  O3' C2' C4' H3'
DC  C3' C5' O4' H4'
DC  C4' O5' H5'1 H5'2
#
DC5  C2' O4' N1 H1'
DC5  C3' C1' H2'1 H2'2
DC5  O3' C2' C4' H3'
DC5  C3' C5' O4' H4'
DC5  C4' O5' H5'1 H5'2
#
DC3  C2' O4' N1 H1'
DC3  C3' C1' H2'1 H2'2
DC3  O3' C2' C4' H3'
DC3  C3' C5' O4' H4'
DC3  C4' O5' H5'1 H5'2
#
DT  C2' O4' N1 H1'
DT  C3' C1' H2'1 H2'2
DT  O3' C2' C4' H3'
DT  C3' C5' O4' H4'
DT  C4' O5' H5'1 H5'2
#
DT5  C2' O4' N1 H1'
DT5  C3' C1' H2'1 H2'2
DT5  O3' C2' C4' H3'
DT5  C3' C5' O4' H4'
DT5  C4' O5' H5'1 H5'2
#
DT3  C2' O4' N1 H1'
DT3  C3' C1' H2'1 H2'2
DT3  O3' C2' C4' H3'
DT3  C3' C5' O4' H4'
DT3  C4' O5' H5'1 H5'2
#
DG  C2' O4' N9 H1'
DG  C3' C1' H2'1 H2'2
DG  O3' C2' C4' H3'
DG  C3' C5' O4' H4'
DG  C4' O5' H5'1 H5'2
#
DG3  C2' O4' N9 H1'
DG3  C3' C1' H2'1 H2'2
DG3  O3' C2' C4' H3'
DG3  C3' C5' O4' H4'
DG3  C4' O5' H5'1 H5'2
#
DG5  C2' O4' N9 H1'
DG5  C3' C1' H2'1 H2'2
DG5  O3' C2' C4' H3'
DG5  C3' C5' O4' H4'
DG5  C4' O5' H5'1 H5'2
#
DA  C2' O4' N9 H1'
DA  C3' C1' H2'1 H2'2
DA  O3' C2' C4' H3'
DA  C3' C5' O4' H4'
DA  C4' O5' H5'1 H5'2
#
DA3  C2' O4' N9 H1'
DA3  C3' C1' H2'1 H2'2
DA3  O3' C2' C4' H3'
DA3  C3' C5' O4' H4'
DA3  C4' O5' H5'1 H5'2
#
DA5  C2' O4' N9 H1'
DA5  C3' C1' H2'1 H2'2
DA5  O3' C2' C4' H3'
DA5  C3' C5' O4' H4'
DA5  C4' O5' H5'1 H5'2
eof
$AWK '
BEGIN {
	while ( getline < "chiral_defs" > 0 ) {
		if (substr($1,1,1) != "#") {
			icons++
			rescons[icons] = $1
			atcons1[icons] = $2
			atcons2[icons] = $3
			atcons3[icons] = $4
			atcons4[icons] = $5
		}
	}
	ncons = icons
}
$1=="ATOM" {
	name = $3; n1 = substr($3,1,1)
	if( n1=="1" || n1=="2" || n1=="3" ){ name = substr($3,2,length($3)-1) n1 }
	ind = name "," $5
	atno[ind] = $2; restype[$5] = $4
	lastres = $5
}
END { OFS = ""
		
	for (ires=1; ires<=lastres; ires++ ) {
		for (icons=1; icons<=ncons; icons++) {
			if (rescons[icons]==restype[ires]) {
				ind1 = atcons1[icons] "," ires	
				ind2 = atcons2[icons] "," ires	
				ind3 = atcons3[icons] "," ires	
				ind4 = atcons4[icons] "," ires	
				if (ind1 in atno && ind2 in atno && ind3 in atno \
					&& ind4 in atno ) {
					print "#"
					print "#  chirality for residue ",ires, " atoms:  ", \
						atcons1[icons]," ",atcons2[icons]," ",atcons3[icons], \
						" ",atcons4[icons]
					icount++
					if (icount==1) {
						print " &rst iat=",atno[ind1],",",atno[ind2],",", \
						atno[ind3],",",atno[ind4],", "
						print "   r1=10., r2=60., r3=80., r4=130.,"\
						" rk2 = 10., rk3=10.,  &end"
					} else {
						print " &rst iat=",atno[ind1],",",atno[ind2],",", \
						atno[ind3],",",atno[ind4],", &end"
					}
				} else {
					print "bad constraint for residue ",ires, " atoms:  ", \
						atcons1[icons]," ",atcons2[icons]," ",atcons3[icons], \
						" ",atcons4[icons] 
				}
			}
		}
	}
	icount = 0
	for (ires=2; ires<=lastres; ires++ ) {
		inda = "CA," ires
		indn = "N," ires; indc = "C," ires
		iresm1 = ires-1
		indcm = "C," iresm1; indam = "CA," iresm1
		if (inda in atno && indn in atno && indcm in atno && indam in atno) {
			icount++
			print "#"
			print "#  trans-omega constraint for residue ",ires
			if (icount==1) {
				print " &rst iat=",atno[inda],",",atno[indn],",", \
					atno[indcm],",",atno[indam],", "
				print "   r1=150., r2=170., r3=190., r4=210.,"\
					" rk2 = 50., rk3=50.,  &end"
			} else {
				print " &rst iat=",atno[inda],",",atno[indn],",", \
					atno[indcm],",",atno[indam],", &end"
			}
		}
	}
}'  $1 > $2
#
#  tip user off about results & clean up tmp file
#
COUNT=`grep rst $2 | wc -l `
echo $COUNT chiral restraints generated in $2
/bin/rm -f chiral_defs