// --------------------------------------------------------------------
// wm3d - A Flash Molecular Viewer
//
// Copyright (c) 2011-2013, tamanegi (tamanegi@users.sourceforge.jp)
// All rights reserved.
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
// --------------------------------------------------------------------

package pdb;

import tinylib.Point3D;

class WMHydrogenBond {
  /**
    arbitrary hydrogen bond energy
  **/
  @:isVar public var E( get, set ):Float;
    /**
      getter of `E`
    **/
    public function get_E():Float { return( E ); }
    /**
      setter of `E`
    **/
    public function set_E( e:Float ):Float {
      E = e;
      return( e );
    }

  /**
    energy component: distance part
  **/
  @:isVar public var Er( get, set ):Float;
    /**
      getter of `Er`
    **/
    public function get_Er():Float { return( Er ); }
    /**
      setter of `Er`
    **/
    public function set_Er( e:Float ):Float {
      Er = e;
      return( Er );
    }
  /**
    energy component: in-plane angle part
  **/
  @:isVar public var Ep( get, set ):Float;
    /**
      getter of `Ep`
    **/
    public function get_Ep():Float { return( Ep ); }
    /**
      setter of `Ep`
    **/
    public function set_Ep( e:Float ):Float {
      Ep = e;
      return( Ep );
    }
  /**
    energy component: out-of-plane angle part
  **/
  @:isVar public var Et( get, set ):Float;
    /**
      getter of `Et`
    **/
    public function get_Et():Float { return( Et ); }
    /**
      setter of `Et`
    **/
    public function set_Et( e:Float ):Float {
      Et = e;
      return( Et );
    }

  // -----------------------------------------------------------------------
  /**
    Constructor.
  **/
  public function new() {
    Er = 0.0;
    Ep = 0.0;
    Et = 0.0;
    E = 0.0;
  }

  /**
    calculate hydrogen bond energy from pre-determined `Er, `Ep`, and `Et`.
  **/
  public function calcE() {
    E = Er * Ep * Et;
  }

  /**
    calculate distance part energy, where `dno` is the length parameter
  **/
  public function calcEr( dno:Float ):Float {
    var d:Float = HBParams.HB_dist / dno;
    if ( d > 1.0 ) {
      Er = HBParams.HB_ene;
    } else {
      // Er = D / r^6 + C / r^8
      var d2:Float = d * d;
      var d6:Float = d2 * d2 * d2;
      var d8:Float = d2 * d6;
      // there may be misprint in the paper
      // C = -3 Em rm^8, D = -4 Em rm^6 in the paper.
      // this maybe C = -3 Em rm^8, D = 4 Em rm^6
      Er = HBParams.HB_ene * ( 4.0 * d6 - 3.0 * d8 );
    }
    return( Er );
  }

  /**
    calculate in-plane angle part of HB energy for X-H-Y.
  **/
  public function calcEp( vho:Point3D,
                          vnh:Point3D,
                          dho:Float,
                          dnh:Float ):Float {
    // Ep = cos^2 0.5p, where p is the angle of N-H-O
    var cosangle:Float = - vho.dot( vnh ) / ( dho * dnh );
    if ( cosangle > 0.0 ) return( 0.0 );
    var angle:Float = Math.PI - Math.acos( cosangle );
    cosangle = Math.cos( 0.5 * angle );
    Ep = cosangle * cosangle;
    return( Ep );
  }

  /**
    calculate out-of-plane angle part of hydrogen bond energy
  **/
  public function calcEt( posH:Point3D,
                          posO:Point3D,
                          posC:Point3D,
                          posCA:Point3D ):Float {
    var vc_ca = Point3D.getSub( posCA, posC );
    var vc_o = Point3D.getSub( posO, posC );
    var vertical = Point3D.getCross( vc_ca, vc_o );
    vertical.normalize();
    var vo_h = Point3D.getSub( posH, posO );
    var vo_h_plane = Point3D.getSub( vo_h, vertical.multiply( vertical.dot( vo_h ) ) );
    var costi = vo_h.dot( vc_o ) / ( vo_h.norm() * vc_o.norm() );
    var ti = Math.acos( costi );
    var costo = vo_h.dot( vo_h_plane ) / ( vo_h.norm() * vo_h_plane.norm() );
    // there maybe misprint; the second Ep = equation should be Et =
    // note: why 45 degrees point is the most stable for ti?
    if ( ti > Math.PI / 3.0 ) {
      // decrease steeply
      //Et = costo * Math.pow( Math.cos( ti - Math.PI / 3.0 ), 4 );
      Et = costo * Math.cos( ti - Math.PI / 3.0 );
    } else {
      // moderate decrease
      //Et = costo * Math.pow( Math.cos( ti - Math.PI / 3.0 ), 2 );
      Et = costo * Math.pow( Math.cos( ti - Math.PI / 3.0 ), 0.5 );
    }
    return( Et );
  }
}

/**
  Container for hydrogen bond energy between two residues.
  This only consider hydrogen bonds between backbone N-H and C=O.
**/
class BackboneHB {
  /**
    hydrogen bond energy for res0(N) - res1(O)
  **/
  public var E01:WMHydrogenBond;
  /**
    hydrogen bond energy for res1(N) - res0(O)
  **/
  public var E10:WMHydrogenBond;

  // ------------------------------------------------------------------------
  /**
    Constructor. If two residues are given, backbone hydrogen bond energy
    is calculated. See also `calc` function.
  **/
  public function new( ?r0:PdbResidue,
                       ?r1:PdbResidue ) {
    E01 = new WMHydrogenBond();
    E10 = new WMHydrogenBond();
    if ( r0 != null && r1 != null ) calc( r0, r1 );
  }

  /**
    calculate backbone hydrogen bond between two residues.
    Note that this function does not return hydrogen bond energy.
    Use `getEne` function if you need hydrogen bond energy or `bonded`
    function if you want to know whether hydrogen bond exists between residues.
  **/
  public function calc( r0:PdbResidue,
                        r1:PdbResidue ) {
    if ( !r0.isAmino() || !r1.isAmino() ) return;
    if ( r0.getAtom( "H" ) == null && r1.getAtom( "H" ) == null ) return;
    if ( r0.getAtom( "H" ) != null ) {
      // 0(NH) -> H(OC)
      var vno = Point3D.getSub( r0.getAtom( "N" ).pos, r1.getAtom( "O" ).pos );
      var distNO = vno.norm();
      var vho = Point3D.getSub( r0.getAtom( "H" ).pos, r1.getAtom( "O" ).pos );
      var distHO = vho.norm();
      var vnh = Point3D.getSub( r0.getAtom( "N" ).pos, r0.getAtom( "H" ).pos );
      var distNH = vnh.norm();
      E01.calcEr( distNO );
      E01.calcEp( vho, vnh, distHO, distNH );
      E01.calcEt( r0.getAtom( "H" ).pos, r1.getAtom( "O" ).pos,
                  r1.getAtom( "C" ).pos, r1.getAtom( "CA" ).pos );
      E01.calcE();
    }
    if ( r1.getAtom( "H" ) != null ) {
      // 1(NH) -> 0(OC)
      var vno = Point3D.getSub( r1.getAtom( "N" ).pos, r0.getAtom( "O" ).pos );
      var distNO = vno.norm();
      var vho = Point3D.getSub( r1.getAtom( "H" ).pos, r0.getAtom( "O" ).pos );
      var distHO = vho.norm();
      var vnh = Point3D.getSub( r1.getAtom( "N" ).pos, r1.getAtom( "H" ).pos );
      var distNH = vnh.norm();
      E10.calcEr( distNO );
      E10.calcEp( vho, vnh, distHO, distNH );
      E10.calcEt( r1.getAtom( "H" ).pos, r0.getAtom( "O" ).pos,
                  r0.getAtom( "C" ).pos, r0.getAtom( "CA" ).pos );
      E10.calcE();
    }
  }

  /**
    return hydrogen bond energy in this definition
  **/
  public function getEne():Float {
    return( Math.min( E01.E, E10.E ) );
  }

  /**
    return whether there is a hydrogen bond or not
  **/
  public function bonded( threshold:Float ):Bool {
    if ( getEne() < threshold ) return( true );
    return( false );
  }
}
