#!/usr/local/bin/gawk -f
#
#   
#@   purpose: create restraints for working with proxy residues
#
#   input: seq
#   
#   

BEGIN { 
    stderr="/dev/stderr"

    sel_duh = "C* / HBOUND - @HIS* @PHE @TRP @TYR CD* CE* CZ* CH* - @*DU*"
    sel_dub = "C* - METHYL / HBOUND - @TYR @PHE @TRP @HIS CD* CE* CZ* CH* - CA + @GLY CA - @*DU*"
    sel_dut = "C* / METHYL / HBOUND - @*DU*"
    sel_dur = "C* / METHYL / HBOUND  - @*DU*"
    sel_duf = "@PHE @TYR @TRP @HIS CD* CE* CZ* CH* / HBOUND"
    sel_dua = "@PHE @TYR @TRP @HIS CD* CE* CZ* CH* / HBOUND"
    sel_dun = "N* / HBOUND - @*DU*"
    sel_xdun= "@ASN ND2 + @GLN NE2  - @*DU*"
    sel_qdux= "C* - @*DU* - @TYR @PHE @TRP @HIS CD* CE* CZ* CH* - @GLY CA "
    sel_xcc = "C* - @*DU* - @TYR @PHE @TRP @HIS CD* CE* CZ* CH* - @GLY CA "
    
    sel_asg_ch1 = "CA + @VAL @THR @ILE CB + @LEU CG / HBOUND - @GLY CA - @*DU* -"
    sel_asg_ch2 = "C* - METHYL -BB / HBOUND + @GLY CA - @VAL @THR @ILE CB - @LEU CG - @TYR @PHE @TRP @HIS CD* CE* CZ* CH* - @*DU* - CA + @GLY CA"
    sel_asg_ch3 = "C* / METHYL / HBOUND - @*DU*"
    sel_asg_cha = "@PHE @TYR @TRP @HIS CD* CE* CZ* CH* / HBOUND"
    sel_asg_nh1 = "N* / HBOUND - @ASN @GLN ND2 NE2 - @*DU*"
    sel_asg_nh2 = "@ASN ND2 + @GLN NE2"
    sel_asg_nhw = "@TRP NE1"
    
}
    #prot_selX="CA CB CG* CD* CE* / HBOUND - @*DU*"
    #prot_selX="CA + @VAL @THR @ILE CB + @LEU CG / HBOUND - @*DU*"
    #prot_selX="C* / HBOUND - @*DU* - @HIS* @PHE @TRP @TYR CD* CE* CZ* CH*"
    #prot_selX="CA + CB - @ALA + CG* - @VAL + CD* - @LEU + CE* - @MET  / HBOUND - @*DU*"
    #prot_selX="C* - METHYL  / HBOUND - @TYR @PHE @TRP @HIS CD* CE* CZ* CH* - @*DU* - CA + @GLY CA"
    #prot_selX="@ALA CB + @VAL CG* + @ILE @THR CG2 + @LEU @ILE CD* + @MET CE / HBOUND - @*DU*"
    #prot_selX="C* / METHYL / HBOUND  - @*DU*"
    #prot_selX="@PHE @TYR CD* CE* CZ* + @TRP CD1 CE* CZ* CH* + @HIS* CD2 CE1 / HBOUND - @*DU*"
    #prot_selX="@PHE @TYR CD* CE* CZ* + @TRP CD1 CE* CZ* CH* + @HIS* CD2 CE1 / HBOUND  - @*DU*"
    #prot_selX="@PHE @TYR CD* CE* CZ* + @TRP CD1 CE* CZ* CH* + @HIS* CD2 CE1 / HBOUND  - @*DU*"
    #prot_selX=" N + @TRP NE NH* + @ARG NE NH* + @LYS NZ + @ASN ND2 + @GLN NE2 + @HIS ND1 NE2 - @*DU*"
    # prot_selX=" @ASN ND2 + @GLN NE2  - @*DU*"

$1~/DUH$/ { pduh++  # 
    dummy_heavyatoms = "CA"
    weight=2.0
    prot_selX=sel_duh
}
$1~/DUB$/ { pdub++  # 
    dummy_heavyatoms = "CA"
    weight=2.0
    prot_selX=sel_dub
}

$1~/DUT$/ { pdut++  # 
    dummy_heavyatoms = "CA"
    weight=2.0
    prot_selX=sel_dut
}

$1~/DUR$/ { pdur++  # 
    dummy_heavyatoms = "CD1 CD2 CE1 CE2 CZ*"
    weight=0.5
    prot_selX=sel_dur
}

$1~/DUF$/ { pduf++  # 
    dummy_heavyatoms = "CD*"  # "CD1 CD2"
    weight=1.0
    prot_selX=sel_duf
}

$1~/DUA$/ { pdua++  # 
    dummy_heavyatoms = "CD*"
    weight=1.0
    prot_selX=sel_dua
}

$1~/DUN$/ { pdun++  # 
    dummy_heavyatoms = "N"
    weight=2.0
    prot_selX=sel_dun
}
$1~/XDUN$/ { pdun++  # 
    dummy_heavyatoms = "N"
    weight=2.0
    
    prot_selX=sel_xdun
}
$1~/DUX$/ { pdux++  # 
    dummy_heavyatoms = "CA CB CG1 CG2 CD1"
    weight=2.0
    
    prot_selX=sel_qdux
}
$1~/^XCC$/ { pxcc++  # 
    dummy_heavyatoms = "CA CB"
    weight=2.0
    
    prot_selX=sel_xcc
}

/DU/{ 
    weight=2.0
    if (/#(asg|select)=/) {
        #print > stderr
        if      (/#asg=CH1/) prot_selX=sel_asg_ch1
        else if (/#asg=CH2/) prot_selX=sel_asg_ch2
        else if (/#asg=CH3/) prot_selX=sel_asg_ch3
        else if (/#asg=CHA/) prot_selX=sel_asg_cha
        else if (/#asg=NH1/) prot_selX=sel_asg_nh1
        else if (/#asg=NH2/) prot_selX=sel_asg_nh2
        else if (/#asg=NHW/) prot_selX=sel_asg_nhw
        else {
            if (/#select=/) prot_selX=gensub(/.*select=["]([^"]+)["].*/,"\\1","")
            print $0,"not recognized" > stderr
        }
    }
    source[$2]=$0
    make_dist(prot_selX,$2,dummy_heavyatoms,"0.01",weight)
    next
}

END {
    #
    # generate lower bounds for repulsion
    #
    lol=1.5
    weight=10.0    
    
    print "\n# lower bound restraints to replace vdW betweeen dummies" 
    #print "distance make \"@*DU* CA CD* CE* CZ* N\"",
    #        " \"@*DU* CA CD* CE* CZ* N\" lol="lol" weight="weight
    print "distance make \"@*DU* CA CD* CE* CZ* N / ASSIGNED\"",
            " \"@*DU* CA CD* CE* CZ* N / ASSIGNED\" lol="lol" weight="weight
    
    print "write upl id"  
    print "write lol du_vdw"
}

function make_dist(protCs,ndumres,dummyCs,dupl,weight,   ndumC,ddumC,d,dmformat) { 
    ndumC=split(dummyCs,ddumC)
    dmformat="distance make first=\"%i %s\" second=\"%s\" upl=%s weight=%s ambiguous\n"
    if (dupl=="")   dupl=0.01
    if (weight=="") weight=2.0
    if (ndumC>1)print ""
    printf "# %s %s\n",ndumres,source[ndumres]
    for(d=1;d<=ndumC;d++) {
        printf dmformat,ndumres,ddumC[d],protCs,dupl,weight/ndumC
    }
    return ndumC

}