#!/usr/bin/perl
# Purpose: read in a PDB NMR restraint dat file and translate the NOE/hydrogen
# bonding and dihedral torsion restraints into CHARMM input files.
#
# 1. NOE restraints
#
# === MR format:
#  ...
#  assign (resid   174 and name HN ) (resid   171 and name HB2 )) 2.439  0.714 0.714
#  ....
# ==== CHARMM format (NEW NOE with soft asypmptote)
#  NOE
#    ...
#    ASSI rexp 3 kmin 3.00 kmax 1.00 fmax 1.00 rmin 99.00 rmax 0.00 rswi 1.0 sexp 1.0 SUMR -
#         SELE (resid   174 .AND. type HN) END -
#         SELE (resid   174 .AND. type HN) .OR. (...) -
#              .OR. (...) END
#  END
#
# NOTE: kmin=1.0 kmax=1.0 can only be changed from the source code.
#       (The right way is to use NOE SCALE command in CHARMM).
# NOTE: rexp and fmax can be specified from command line.
#       If omitted, default of rexp =-1/6, fmax = 99.0 will be used.
#
# 3.HYDROGEN BONDING (implement as NOE distance restraint)
# 
use strict;

sub usage {
    print STDERR 
	"usage: mr2charmm.pl MRFILE CHARMMFILE [options]\n";
    exit(1);
}

# read in command line options
&usage if($#ARGV<1);
my $MRFILE = $ARGV[0];
my $CHARMMFILE = $ARGV[1];
shift @ARGV;
shift @ARGV;

my $kmin = 1.0;         # NOE restraints parameters
my $kmax = 1.0;
my $fmax = 1.0;
my $rswi = 3.0;
my $sexp = 1.0;
my $rexp = -1.0/6;          # MR default averaging for ambiguous NOEs
while ($#ARGV >=0 ) {   # optional fmax, rexp specifications
    if ( $ARGV[0] =~ /fmax=([+-]*[\d\.]+)/ ) {
	$fmax = $1;
	shift @ARGV;
    } elsif ( $ARGV[0] =~ /rexp=([+-]*[\d\.]+)/ ) {
	$rexp = $1;
	shift @ARGV;
    } else {
	print STDERR "\n>>> Invalid option: $ARGV[0]. Skipped.\n";
	shift @ARGV;
    }
}

my $noeparam = "
if \@\?rexp eq 0   set rexp = $rexp
if \@\?kmin eq 0   set kmin = $kmin
if \@\?kmax eq 0   set kmax = $kmax
if \@\?fmax eq 0   set fmax = $fmax
if \@\?rswi eq 0   set rswi = $rswi
if \@\?sexp eq 0   set sexp = $sexp
";

# open files for input/output
open(MR, "<$MRFILE") || 
    die "Error: can not open file $MRFILE for read.\n";
print("Open an MR file: $MRFILE for read.\n");

my $AllNOEFile = "$CHARMMFILE.noe.all";             # all
my $IntraNOEFile = "$CHARMMFILE.noe.intra";         # |i-j|=0
my $SeqNOEFile = "$CHARMMFILE.noe.seq";             # |i-j|=1
my $MediumNOEFile = "$CHARMMFILE.noe.med";          # |i-j|<=4
my $LongNOEFile = "$CHARMMFILE.noe.long";           # |i-j|>4
open(AllNOE, ">$AllNOEFile") || 
    die "Error: can not open file $AllNOEFile for write.\n";
open(IntNOE, ">$IntraNOEFile") || 
    die "Error: can not open file $IntraNOEFile for write.\n";
open(SeqNOE, ">$SeqNOEFile" ) || 
    die "Error: can not open file $SeqNOEFile for write.\n";
open(MedNOE, ">$MediumNOEFile") || 
    die "Error: can not open file $MediumNOEFile for write.\n";
open(LngNOE, ">$LongNOEFile") || 
    die "Error: can not open file $LongNOEFile for write.\n";
print "
NOE restraints are grouped and written to CHARMM input files:

                   All NOE: $AllNOEFile
          Intraresidue NOE: $IntraNOEFile
             Sequetial NOE: $SeqNOEFile
          Medium range NOE: $MediumNOEFile
            Long range NOE: $LongNOEFile
";
print ("\nNOE restraints parameters: \n");
print ("   kmin = $kmin, kmax = $kmax, rexp = $rexp, 
   fmax = $fmax, rswi = $rswi, sexp = $sexp, SUMR is on\n\n");

# writing headers
my $pwd = `pwd`;
chop($pwd);
my $date = `date`;
chop($date);
my $user = `whoami`;
chop($user);
my $header = "* From NMR restraint file: $pwd/$MRFILE
* created on: $date \t by user: $user\n";

printf AllNOE $header;
printf AllNOE "* All NOEs \n*\n";
printf AllNOE $noeparam;
printf AllNOE "NOE\n";

printf IntNOE $header;
printf IntNOE "* Intraresidue NOEs \n*\n";
printf IntNOE $noeparam;
printf IntNOE "NOE\n";

printf SeqNOE $header;
printf SeqNOE "* Sequential NOEs \n*\n";
printf SeqNOE $noeparam;
printf SeqNOE "NOE\n";

printf MedNOE $header;
printf MedNOE "* Medium range NOEs \n*\n";
printf MedNOE $noeparam;
printf MedNOE "NOE\n";

printf LngNOE $header;
printf LngNOE "* Long range NOEs \n*\n";
printf LngNOE $noeparam;
printf LngNOE "NOE\n";

# start reading MRFILE
my $newline; 
my $blocktype = "noe";
my($nintnoe,$nseqnoe,$nmednoe,$nlngnoe,$ndihe);
$nintnoe=0,$nseqnoe=0,$nmednoe=0,$nlngnoe=0,$ndihe=0;
while(<MR>) {
    $newline = $_; 
    &readNOE;
}

my $ntotalnoe = $nintnoe + $nseqnoe + $nmednoe + $nlngnoe; 
printf "
------------- Conversion summary ----------------- 

         total number of NOEs: $ntotalnoe

            intraresidue NOEs: $nintnoe
              sequential NOEs: $nseqnoe
            medium range NOEs: $nmednoe
              long range NOEs: $nlngnoe

\n";
printf AllNOE "END\nRETURN\n! total number of restraints: $ntotalnoe\n";
printf IntNOE "END\nRETURN\n! total number of restraints: $nintnoe\n";
printf SeqNOE "END\nRETURN\n! total number of restraints: $nseqnoe\n";
printf MedNOE "END\nRETURN\n! total number of restraints: $nmednoe\n";
printf LngNOE "END\nRETURN\n! total number of restraints: $nlngnoe\n";

close(MR);
close(IntNOE);
close(SeqNOE);
close(MedNOE);
close(LngNOE);
exit(0);

#
# read in NOE block, group the assignments according to |i-j|
#   atom selection format is defined by $atomsel (has to be a single line)
# ...
# assign (resid  80 and name HB#  ) (resid  80 and name HN  )  3.1  1.1  1.1
# assign (resid  82 and (name HB1 or name HB2  )) (resid  82 and name HN  )  2.8  1.0  1.0
# assign (resid  85 and (name HE21 or name HE22  )) (resid  85 and (name HG1 or name HG2  )) 2.9  1.0  1.0
# assign (resid  85 and name HB#  ) (resid  85 and (name HE21 or name HE22  )) 3.5  1.2  1.2
# ...
#
sub readNOE {
    my $atomsel = "\\(\\s*resi(d|)\\s+(\\d+)\\s+and\\s+([\\s\\S]+)\\s*\\)";
    my($resid1,$resid2,$atom1,$atom2,$rmean,$rplus,$rminus,$rmin,$rmax);
    my $assi;

    if($newline=~/^\s*assi[a-z]*/) {  # start of an assignment	
	if($newline =~ /$atomsel\s*$atomsel\s*([\d\.]+)\s+([\d\.]+)\s+([\d\.]+)/i) { #single line
	    $resid1 = $2, $atom1=$3, $resid2 = $5, $atom2 = $6;
	    $rmean = $7, $rminus = $8, $rplus = $9;	    

	    $atom1 =~ s/or/.or./ig;
	    $atom1 =~ s/name/type/ig;
	    $atom2 =~ s/or/.or./ig;
	    $atom2 =~ s/name/type/ig;
	    $rmin = $rmean - $rminus;
	    $rmax = $rmean + $rplus;
	    $assi="  ASSI rmin $rmin rmax $rmax -
       SELE resid $resid1 .AND. $atom1 end SELE resid $resid2 .AND. $atom2 end -
       rexp \@rexp fmax \@fmax rswi \@rswi sexp \@sexp kmin \@kmin kmax \@kmax SUMR\n";

	    if(abs($resid1-$resid2)==0) {
		printf IntNOE $assi;
		$nintnoe++;
	    } elsif(abs($resid1-$resid2)==1) {
		printf SeqNOE $assi;
		$nseqnoe++;
	    } elsif(abs($resid1-$resid2)<=4) {
		printf MedNOE $assi;
		$nmednoe++;
	    } else {
		printf LngNOE $assi;
		$nlngnoe++;
	    }
	    printf AllNOE $assi;

#	} elsif ($newline =~ /$atomsel/i) {  #double line
	    
	} else {
	    printf  ">>> $newline Error: unrecogonized format\n";
	}	
    } else {	
	print IntNOE "! $newline";
	print AllNOE "! $newline";
    }
    
    return;
}
