#!/usr/local/bin/perl -w
=pod

dihedrals2xplor v. 14 Jun 1996, Dave Schweisguth <dcs@proton.chem.yale.edu>

Converts a table of dihedral angles to Xplor dihedral restraints format.
Idea from Jon Lapham's cdih_make.

Input file format:

The input file is a series of lines, each beginning with
- a residue (a, c, g, t or u, possibly followed by other letters, which will be
  ignored), optionally followed by a residue number (no whitespace between
  nucleotide and number), followed by
  - nothing, to restrain each angle in that residue to the current default
  - "none", to leave that residue entirely unconstrained
  - eleven angle restraints, in the order alpha, beta, gamma, epsilon, zeta,
    chi, nu0-nu4, to restrain that residue to those values. Each restraint may
    be a pair of numbers, the first being the value and the second the range,
    or "none", to leave that individual angle entirely unconstrained. Specify
    eleven restraints even for a nucleotide which does not have eleven dihedral
    angles (e.g. 5' alpha or 3' epsilon or zeta); use "none" for nonexistent
    angles.
- "s[egment]" and a segid
- "r[esidue]" and a residue number
- "d[efault]", followed by
  - nothing, to reset the defaults to the original defaults included in the
    program
  - eleven angle restraints, just as for a residue, to be used for subsequent
    residues which do not specify individual restraints

For example,

---snip---
# ITALY2 stem, all A-form defaults
segment a
g
g
c
a
g
g
g
segment b
c
c
c
u
g
c
c
---snip---

or

---snip---
# Dihedral angle constraints for ITALY2
# A-form values from Arnott and Chandrasekaran (1982) as reported by Sanger (1984)

#   alpha    beta    gamma  epsilon   zeta    chi      nu0       nu1        nu2      nu3       nu4

g1  none     none    none   -153  15  -71 15  -158 15  -27.5 15    7.7 15   15.0 15  -32.0 15  36.8 15
g2   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c3   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
a4   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
g5   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
g6   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
g7   -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c8   -68 30  178 30  none   -153 120  -71 30  -158 30    3   30  -25   30   37   30  -36   30  21   30
u9   -68 30  178 30  none   -153 120  none    -158 30    3   30  -25   30   37   30  -36   30  21   30
c10  none    178 30  54 30  -153 120  none    -158 30  -36.0 30   24.2 30  -3.2  30  -19.1 30  34.1 30
a11  none    178 30  none   -153 120  -71 30  -158 30  -34.3 30   19.6 30   2.6  30  -23.8 30  35.9 30
u12  -68 30  178 30  none   -153 120  none    -158 30  -35.2 30   35.2 30  -21.7 30    0.0 30  21.8 30
a13  none    178 30  none   -153 120  -71 30  -158 30  -29.4 30   10.5 30   12.3 30  -30.5 30  37.0 30
a14  -68 30  178 30  54 30  -153 120  -71 30  -158 30    3   30  -25   30   37   30  -36   30  21   30
c15  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c16  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c17  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
u18  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
g19  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c20  -68 15  178 15  54 15  -153  15  -71 15  -158 15    3   15  -25   15   37   15  -36   15  21   15
c21  -68 15  178 15  54 15  none      none    -158 15  -32.5 15   15.9 15    6.8 15  -26.9 15  36.7 15
---snip---

or (the same thing, using defaults)

---snip---
# Dihedral angle constraints for ITALY2
# A-form values from Arnott and Chandrasekaran (1982) as reported by Sanger (1984)

#  alpha   beta    gamma  epsilon   zeta    chi      nu0       nu1        nu2      nu3       nu4

g  none    none    none   -153  15  -71 15  -158 15  -27.5 15    7.7 15   15.0 15  -32.0 15  36.8 15
g
c
a
g
g
g
# Increase range to +/- 30 degrees in the loop
d  -68 30  178 30  54 30  -153 120  -71 30  -158 30    3   30  -25   30   37   30  -36   30  21   30
c  -68 30  178 30  none   -153 120  -71 30  -158 30    3   30  -25   30   37   30  -36   30  21   30
u  -68 30  178 30  none   -153 120  none    -158 30    3   30  -25   30   37   30  -36   30  21   30
c  none    178 30  54 30  -153 120  none    -158 30  -36.0 30   24.2 30  -3.2  30  -19.1 30  34.1 30
a  none    178 30  none   -153 120  -71 30  -158 30  -34.3 30   19.6 30   2.6  30  -23.8 30  35.9 30
u  -68 30  178 30  none   -153 120  none    -158 30  -35.2 30   35.2 30  -21.7 30    0.0 30  21.8 30
a  none    178 30  none   -153 120  -71 30  -158 30  -29.4 30   10.5 30   12.3 30  -30.5 30  37.0 30
a
# Reduce range back to +/- 15 degrees in the stem
d
c
c
c
u
g
c
c  -68 15  178 15  54 15  none      none    -158 15  -32.5 15   15.9 15    6.8 15  -26.9 15  36.7 15
---snip---

Blank lines and lines beginning with # or ; are ignored. Lines beginning with
! are printed, thus becoming Xplor comments in the output file.

Segment name definitions are optional. Without them, the entire molecule is
segment a. Defining a new segment name resets the residue number to 1.

Residue number definitions are optional. They may be negative. Without them,
numbering begins at 1.

Caveats:
- Assumes residue numbering to be continuous within each segment
- Assumes that a nucleotide at the 5' end of a segment has a 5' phosphate
  (which does not need an alpha constraint, but does need a beta constraint)
  and a nucleotide at the 3' end of a segment has a 3' hydroxyl (which does not
  need epsilon or zeta constraints). If you want to define the equivalent of
  alpha for a 5' nucleotide or epsilon or zeta for a 3' nucleotide, you'll have
  to do it by hand. If your 5' nucleotide has a 5' hydroxyl, specify its beta
  angle to be "none" or remove the definition from the output by hand.

=cut

### Preliminaries

require 5.002;			    # Perl 5.002 required
use strict;			    # Require optional-but-desirable practices

### Parameters

# Environment

(my $whatami = $0) =~ s|.*/||;	    # `basename $0`

# Configuration

my $seg = 'a';			    # Default segment name
my $res = 1;			    # Default first residue number
my $segments = 0;		    # Include segment in atom selection?
				    # -1 = never, 0 = maybe, 1 = always
my $types = 'acgtu';		    # Legal nucleotides

my @angle_names = qw(alpha beta gamma epsilon zeta chi nu0 nu1 nu2 nu3 nu4);
				    # Angle names in desired order

# Atoms for each angle

my %atoms;
%{$atoms{a}} = (
    alpha   => ["O3'", "P  ", "O5'", "C5'"],
    beta    => ["P  ", "O5'", "C5'", "C4'"],
    gamma   => ["O5'", "C5'", "C4'", "C3'"],
    epsilon => ["C4'", "C3'", "O3'", "P  "],
    zeta    => ["C3'", "O3'", "P  ", "O5'"],
    chi	    => ["O4'", "C1'", "N9 ", "C4 "],
    nu0	    => ["C4'", "O4'", "C1'", "C2'"],
    nu1	    => ["O4'", "C1'", "C2'", "C3'"],
    nu2	    => ["C1'", "C2'", "C3'", "C4'"],
    nu3	    => ["C2'", "C3'", "C4'", "O4'"],
    nu4	    => ["C3'", "C4'", "O4'", "C1'"]
);
$atoms{g} = $atoms{a};
$atoms{c} = {%{$atoms{a}}};
$atoms{c}{chi} = ["O4'", "C1'", "N1 ", "C2 "];
$atoms{t} = $atoms{u} = $atoms{c};

# A-form values for each angle, from Arnott and Chandrasekaran (1982)
# as reported by Sanger (1984)

my $original_default_angle_values = {
    alpha   => [ -68, 15],
    beta    => [ 178, 15],
    gamma   => [  54, 15],
    epsilon => [-153, 15],
    zeta    => [ -71, 15],
    chi	    => [-158, 15],
    nu0	    => [   3, 15],
    nu1	    => [ -25, 15],
    nu2	    => [  37, 15],
    nu3	    => [ -36, 15],
    nu4	    => [  21, 15],
};

# Initialization (don't change these)

my $segd = 1;
my $segsel = '';
my $mid_segment = 0;
my $default_angle_values = $original_default_angle_values;
my(@input, $line, @line, $angle_name, @resids, $angle_values, $mid_molecule);

### Arguments and error-checking

# Parse args

my($arg, $sign, $first, $rest);
while (@ARGV and ($sign, $first, $rest) = ($ARGV[0] =~ /^([\-+])(.)(.*)/)) {
    if ($sign eq '+' && $first !~ /[s]/) { # -/+ switches (none at the moment)
	&usage("$sign$first is not an option.\n");
    }
    if ($first =~ /[\0]/) {	# Switches with arguments (none at the moment)
    	shift;
    	$arg = $rest ne '' ? $rest : @ARGV ? shift :
      	    &usage("$sign$first requires an argument.\n");
    } elsif ($rest eq '') {
	shift;
    } else {
    	$ARGV[0] = "$sign$rest";
    }
    if    ($first eq 's') { $segments = $sign eq '-' ? 1 : -1; }
    elsif ($first eq 'u') { &usage(0); }
    else		  { &usage("$sign$first is not an option.\n"); }
}

sub usage {
    warn $_[0] ? "$whatami: $_[0]" : '', <<EOP;
Usage: $whatami [-u] [-/+s] [file ...]
-s  Include segment in atom selection	    (default for segmented molecule)
+s  Don't include segment in atom selection (default for unsegmented molecule)
-u  This message
EOP
    exit !! $_[0];
}

# Postprocess arguments

$segsel = "segid $seg and " if $segments == 1;

### Do it

@input = <>;
$line = 0;
foreach (@input) {
    next if /^$/;		    # Skip blank line
    next if /^[#;]/;		    # Skip # or ; comment
    if (/^!/) {			    # Print ! comment
	print;
	next;
    }
    @line = split(' ', $_);
    $first = shift(@line);
    if ($first =~ /^s/i) {	    # segment
    
	shift;
	@line == 1 ||
	    die "$whatami: At $ARGV line $., $first is not followed by a single segid.\n";
	$segd && $mid_molecule &&
	    warn "$whatami: At $ARGV line $., the first segment is defined in mid-molecule.\n";
	$seg = shift @line;
	$segsel = "segid $seg and " unless $segments == -1;
	$segd = 0;
	$res = 1;
	$mid_segment = 0;
	
    } elsif ($first =~ /^r/i) {	    # residue
    
	@line == 1 ||
	    die "$whatami: At $ARGV line $., $first is not followed by a single residue number.\n";
	$res = shift @line;
	&is_int($res) ||
	    die "$whatami: At $ARGV line $., $res is not a number!\n";
	    
    } elsif ($first =~ /^d/i) {	    # default
    
	if (@line == 0) {
	
	    # If 'default' is by itself on the line, use defaults
	
	    $default_angle_values = $original_default_angle_values;
	} else {
	    
	    # The rest of the line must be angle values (or else)
	
	    $default_angle_values = &parse_angle_values;
	}
	
    } elsif ($first =~ /^([$types])(-?\d+)?$/i) {   # acgtu

	if ($2) {		    # Residue number is present
	    $first = $1;
	    $res = $2;
	}

	if (@line == 0) {
	
	    # If the residue is by itself on the line, use defaults
	
	    $angle_values = $default_angle_values;
	} elsif (@line == 1 && $line[0] eq 'none') {
	
	    # If the residue is followed by 'none', don't constrain it
	
	    print "! $first$res is completely unconstrained\n\n";
	    next;
	} else {
	
	    # The rest of the line must be angle values (or else)
	
	    $angle_values = &parse_angle_values;
	}
	
	foreach $angle_name (@angle_names) {

	    if (! $mid_segment && $angle_name eq 'alpha') {
	    
		# Don't constrain alpha in the first residue in a segment
		# (and don't even mention it)

	    } elsif ($angle_name =~ /^epsilon|zeta$/ && &end_segment) {

		# Don't constrain epsilon or zeta in the last residue in a segment
		# (and don't even mention it)

	    } elsif ($$angle_values{$angle_name} eq 'none') {
	    
		# Don't constrain angles set to 'none' (but say so)
	    
		print "! $first$res $angle_name is unconstrained\n\n";
	    } else {

		# Determine the resid for each of the four atoms in this angle
    
		@resids = (
		    $angle_name eq 'alpha' ? $res - 1 : $res,
		    $res,
		    $angle_name eq 'zeta' ? $res + 1 : $res,
		    $angle_name =~ /^epsilon|zeta$/ ? $res + 1 : $res
		);
    
		# Print it
    
		print <<EOP;
! $first$res $angle_name
assign (${segsel}resid $resids[0] and name $atoms{$first}{$angle_name}[0]) (${segsel}resid $resids[1] and name $atoms{$first}{$angle_name}[1])
       (${segsel}resid $resids[2] and name $atoms{$first}{$angle_name}[2]) (${segsel}resid $resids[3] and name $atoms{$first}{$angle_name}[3])
       1 @{$$angle_values{$angle_name}} 2

EOP
	    }
	}

	# Note for later that we've already done at least one residue ...

	$mid_molecule = 1;	    # ... in this molecule
	$mid_segment = 1;	    # ... in this segment
	$res++;
	
    } else {
	die "$whatami: At $ARGV line $., $first means nothing to me.\n";
    }
} continue {
    $line++;			    # Keep track of current line number
}

exit;

### Subroutines

# Return 1 if all args are integers, 0 otherwise

sub is_int {

    foreach (@_) {
	unless (/^-?\d+$/o) {
	    return 0;
	}
    }
    1;

}

# Return 1 if all args are real, 0 otherwise

sub is_real {

    foreach (@_) {
	unless (/^[+-]?(\d*(?:\.\d*)?)(?:e[+-]?\d+)?$/i && $1 ne '.') {
	    return 0;
	}
    }
    1;

}

# Return 1 if at end of segment, 0 otherwise
# Depends entirely on globals ... eeeew.

sub end_segment {
    foreach (@input[$line + 1 .. $#input]) {
	/^[$types]/i && return 0;   # Subsequent line is another residue
	/^s/i && return 1;	    # Subsequent line is a new segment
    }
    1;				    # EOF
}

# Return a pointer to an angle_values structure
# Depends mostly on globals ... eeeew.

sub parse_angle_values {
    my $angle_values = {};

    foreach $angle_name (0 .. $#angle_names) {
    
	@line || die "$whatami: At $ARGV line $., ran out of angle definitions!\n";
    
	if ($line[0] eq 'none') {
	    $$angle_values{@angle_names[$angle_name]} = 'none';
	    shift @line;
	    next;
	}
	
	&is_real($line[0]) ||
	    die "$whatami: At $ARGV line $., $line[0] is not a number!\n";
	&is_real($line[1]) ||
	    die "$whatami: At $ARGV line $., $line[1] is not a number!\n";
	$$angle_values{@angle_names[$angle_name]} = [shift @line, shift @line];
    }
    $angle_values;
}