# -*- coding: utf-8 -*-
# Copyright (c) 2015, PyRETIS Development Team.
# Distributed under the LGPLv3 License. See LICENSE for more info.
"""This file defines the order parameter used for the hydrate example."""
import logging
import os 
import subprocess
import numpy as np
import mdtraj as md
from pyretis.orderparameter import OrderParameter
logger = logging.getLogger(__name__)  # pylint: disable=C0103
logger.addHandler(logging.NullHandler())


class DeadAng1(OrderParameter):
    """RingDiffusion(OrderParameter).

    This class defines the ring diffusion order parameter for
    the hydrate example.

    Attributes
    ----------
    name : string
        A human readable name for the order parameter
    """

    def __init__(self):
        """Initialize the order parameter.

        Parameters
        ----------
        index : tuple of ints
            This is the indices of the atom we will use the position of.
        """
        super().__init__(description='Dihdral angles')
        self.idx1 = np.array([[ 811, 814, 817, 820]
                             ,[ 809, 811, 814, 817]
                             ,[ 807, 809, 811, 814]
                             ,[1655,1657,1659,1662]
                             ,[1657,1659,1662,1663]
                             ,[ 938, 940, 942, 945]
                             ,[ 940, 942, 945, 948]
                             ,[ 942, 945, 948, 949]
                             ,[ 910, 912, 914, 917]
                             ,[ 912, 914, 917, 920]
                             ,[1475,1477,1479,1482]], dtype=np.int32)-1

        self.idx2 = np.array([  820, 821, 822
                              , 823, 824, 825, 826
                              , 827, 828           ], dtype=np.int32)-1


    def calculate(self, system):
        """Calculate the order parameter.

        Here, the order parameter is just the distance between two
        particles.

        Parameters
        ----------
        system : object like :py:class:`.System`
            This object is used for the actual calculation, typically
            only `system.particles.pos` and/or `system.particles.vel`
            will be used. In some cases `system.forcefield` can also be
            used to include specific energies for the order parameter.

        Returns
        -------
        out : float
            The order parameter list.
        """
        file_g96 = system.particles.config[0]
        file_gro = '{}.gro'.format(file_g96)
        cmd = ['gmx', 'editconf', '-f', file_g96, '-o', file_gro]
        exe = subprocess.Popen(cmd, stdin=subprocess.PIPE,
                                   stdout=subprocess.PIPE,
                                   stderr=subprocess.PIPE,shell=False)
        exe.wait()
        return_code = exe.returncode

        trj=md.load(file_gro)
        sasa = md.shrake_rupley(trj)
        total_sasa=0
        for i,area in enumerate(sasa[0]):
            if i in self.idx2: 
                total_sasa += area
  
        orderpALL = md.compute_dihedrals(trj,self.idx1)*180/3.14159
 
        orderpMain =(-orderpALL[0,0]+180)/260
        if orderpMain > 1.1:
            orderpMain = 0
     
        orderp = [orderpMain]
        orderp.extend(orderpALL[0])
        orderp.append(total_sasa)
        print(orderp[0],orderp[1])
        return orderp


class DeadAng2(OrderParameter):
    """RingDiffusion(OrderParameter).

    This class defines the ring diffusion order parameter for
    the hydrate example.

    Attributes
    ----------
    name : string
        A human readable name for the order parameter
    """

    def __init__(self):
        """Initialize the order parameter.

        Parameters
        ----------
        index : tuple of ints
            This is the indices of the atom we will use the position of.
        """
        super().__init__(description='Dihdral angles')
        self.idx1 = np.array([[ 811, 814, 817, 820]
                             ,[ 809, 811, 814, 817]
                             ,[ 807, 809, 811, 814]
                             ,[1655,1657,1659,1662]
                             ,[1657,1659,1662,1663]
                             ,[ 938, 940, 942, 945]
                             ,[ 940, 942, 945, 948]
                             ,[ 942, 945, 948, 949]
                             ,[ 910, 912, 914, 917]
                             ,[ 912, 914, 917, 920]
                             ,[1475,1477,1479,1482]], dtype=np.int32)-1

        self.idx2 = np.array([  820, 821, 822
                              , 823, 824, 825, 826
                              , 827, 828           ], dtype=np.int32)-1


    def calculate(self, system):
        """Calculate the order parameter.

        Here, the order parameter is just the distance between two
        particles.

        Parameters
        ----------
        system : object like :py:class:`.System`
            This object is used for the actual calculation, typically
            only `system.particles.pos` and/or `system.particles.vel`
            will be used. In some cases `system.forcefield` can also be
            used to include specific energies for the order parameter.

        Returns
        -------
        out : float
            The order parameter list.
        """
        file_g96 = system.particles.config[0]
        file_gro = '{}.gro'.format(file_g96)
        cmd = ['gmx', 'editconf', '-f', file_g96, '-o', file_gro]
        exe = subprocess.Popen(cmd, stdin=subprocess.PIPE,
                                   stdout=subprocess.PIPE,
                                   stderr=subprocess.PIPE,shell=False)
        exe.wait()
        return_code = exe.returncode

        trj=md.load(file_gro)
        sasa = md.shrake_rupley(trj)
        total_sasa=0
        for i,area in enumerate(sasa[0]):
            if i in self.idx2: 
                total_sasa += area
  
        orderpALL=md.compute_dihedrals(trj,self.idx1)*180/3.14159
 
        orderpMain = (orderpALL[0,0]+180)/100
        if orderpMain > 1.1:
            orderpMain = 0

        orderp = [orderpMain]
        orderp.extend(orderpALL[0])
        orderp.append(total_sasa)
        return orderp

if __name__ == '__main__':
    testo = DeadAng()
    print('Idx1', testo.idx1[0])