#!/usr/local/bin/gawk -f
#
#@  opposite of nmr2gmx.awk
#@  convert gmx topology distance and orientation restraints 
#
#
#to run do one of: 
# gawk -f gmx2nmr.awk args files......
# gawk -f `which gmx2nmr.awk` args files......
# gmx2nmr.awk args files..... 
# gmx2nmr.awk args files..... 
#
#
#
# -v resoffset=i [default=0]   
#                 offset to match residue number in nmr data file and 
#                 gromacs topology 
#                 nmr-res-nr = gmx-res-nr + resoffset  
#              SO:  resoffset is (pdb) res number at which gmx structure starts -1
#

BEGIN {

	if (!stderr) stderr="/dev/stderr"
    progname="gmx2nmr.awk"
    if ( ARGV[1]~/^(h|help)$/ || help!="" ) { 
        print ARGV[0]
        # because of the which command, help is only available if 
        # this script is in the PATH
        system("gawk '/^#/{print}/^BEGIN/{exit}' `which "progname"`")
        help=1
        exit
    }
    
	mult=0
	if (!functype) functype=1
	if (!weight) weight=1
	if (!type) type="upl"
	# distance larger than maxup nm. are considered lowerbounds ( nonnoes )
	if (!maxup) maxup=2
    if (!default_dico_err) default_dico_err=0.1
    if (!dipco_viol_est) dipco_viol_est=1/2
    
    if (!out) out = "/dev/stdout"
        
    top="top";itp="itp";gro="gro"
}

# get .top of .gro file
FNR==1 {
	if (FILENAME~/\.top$/) { filetype=top; print "#top:",topfile=FILENAME > stderr}
	if (FILENAME~/\.gro$/) { filetype=gro; print "#gro:",grofile=FILENAME > stderr}
	if (FILENAME~/\.itp$/) { filetype=itp; print "#itp:",itpfile=FILENAME > stderr}
}

/[[] *distance[\-_]restraints *[]]/ {
    sub(/[[] +/,"[ ");sub(/ +[]]/," ]")
	filetype=itp
    print "# disre :",disre_file=FILENAME > stderr
}

/[[] *orientation[\-_]restraints *[]]/ {
    sub(/[[] +/,"[ ");sub(/ +[]]/," ]")
	filetype=itp
    print "# orientation restraints :",orire_file=FILENAME > stderr
}

filetype==itp && /; *#/ { # print "###",$0
#/; +dip/ { print "###",$0
    comment=$0;
    if (/^; #/) { sub(/^; #/,"#",comment) ; print comment > out }
    else if (/ ; #/)
        sub(/ ; *#.*/,"",comment)
    # part of comments come originally from dip files; remove that part.
} 

FILENAME==topfile {
	if (/[[] .* \]/) {  itype=$2 ; print $0,itype > stderr }
	if (itype=="atoms") {
		anr=$1
		atom[anr]  = tolower($5)
		resnam[anr]= tolower($4)
		resnum[anr]= tolower($3)
		
		if ( atom[anr]=="h1" && resnum[anr]==1 ) { 
            print $0 > stderr
            atom[anr]="h"
        }
	}

	# adjust resnum for offset between pdb,nmr-data and gmx
	resnum[anr]+= resoffset
}

FILENAME==grofile {
	if (FNR==1) { print "#"$0 > stderr ; next }
	if (FNR==2) { print "#"$0 > stderr ; next }
	if (/SOL/) next

	groform="%5d%5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f" 

	anr = substr($0,16,5)+0
	atom[anr]  = tolower(substr($0,11,5));gsub(/[ \t]/,"",  atom[anr])
	resnam[anr]= tolower(substr($0, 6,5));gsub(/[ \t]/,"",resnam[anr])
	resnum[anr]= substr($0,1,5)+0

    # adjust "h1" in first residue
    if ( atom[anr]=="h1" && resnum[anr]==1 ) atom[anr]=="h"
	#print atomnr" : ",resnum[atomnr],resnam[atomnr],atom[atomnr]

	# adjust resnum for offset between pdb,nmr-data and gmx
	resnum[anr]+= resoffset

}


FILENAME==disre_file || FILENAME==orire_file {
    #print "reading itp-file: ",FILENAME,$0
	if ( /[[] .* \]/ ) { itype=$2 ; print $0,itype > stderr }
	if (/^;/) next
	if (itype=="distance_restraints") {
		if ($1+0<=0) next
		$1=$1
		ai=$1;aj=$2;ftype=$3;rindex=$4;rtype=$5
		low=$6;up1=$7;up2=$8;weight=$9
		amb[rindex]++
		NOE[rindex,amb[rindex]]=ai" "aj" "up1" "low
		if (!atom[ai]) next
		if ( out ) {
			printf  "noerestraint %i %s %s %i %s %s %7.3f %7.3f %i # %s\n",
				resnum[ai],resnam[ai],atom[ai],
				resnum[aj],resnam[aj],atom[aj],
				up1*10,low*10,rindex,$0 > out
		}
	}
	if (itype=="orientation_restraints") {
		if ($1+0<=0) next
		$1=$1
		ai=$1;aj=$2;
        ftype=$3;
        experiment_nr=$4;
        restr_id=$5
		power=$6;
        const=$7;
        dico=$8;
        dico_weight=$9
        
        ## estimate error in dico from weight, this assumes that the weight was
        ## calculated in the same way as in  nmr2gmx.awk
        ## weight = v^2/(v^2+s^2)
        ##     v = estimated violation in Hz 
        ##     s = estimated error in measured  coupling value 
        v2=dipco_viol_est^2   ## = estimated average violation , def=0.5 in nmr2gmx.awk
        if (dico_weight^2 <= 1 ) { # print "#",v2,dico_weight
            s2 = v2/dico_weight-v2
            dico_err = sqrt(s2)
        } else { 
            print  "# don't know how to calculate error in coupling value,"
            printf "# taking %s\n",default_dico_err
            dico_err = default_dico_err
        }
        
        if (!atom[ai]) next
        if (nocomment) comment="" ; else comment=" # "comment
		if ( out ) {
			printf  "dip %4i %4s  %-5s %4i %4s %-5s %8.3f %8.3f %5i%s\n",
				resnum[ai],resnam[ai],atom[ai],
				resnum[aj],resnam[aj],atom[aj],
				dico,dico_err,restr_id,comment > out
		}
	}
}

END {
	if (help)exit
    if (!grofile && !topfile ) {
		print "# no atom indices found , give top or gro file" > stderr
		exit
	}

	imin=imax=rindex
	for ( ri in amb ) { 
		imin = ri+0 <imin ? ri+0 : imin 
	    imax = ri+0 >imax ? ri+0 : imax 
		itot++
	}
    
    if (disre_file!="") {
	    print "# ",itot,"constraints from ",imin,"to",imax 
	    for ( ri=imin;ri<=imax;ri++ ) {
		    for ( p=amb[ri];p>=1;p-- ) {
			    #print NOE[ri,p],p,ri,amb[ri]
			    split(NOE[ri,p],ii)
			    ai=ii[1];aj=ii[2];up1=ii[3];low=ii[4]
			    printf  "constraint %5i %-5s %5s   %5i %-5s %5s   %8.3f %6.3f %3i %5i\n",
				    resnum[ai],tolower(resnam[ai]),tolower(atom[ai]),
				    resnum[aj],tolower(resnam[aj]),tolower(atom[aj]),
				    up1*10,low*10,p,ri 
		    }
	    }
    }
    if (orire_file!="") {  ;
    }
}