# -*- 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 mdtraj as md
import numpy as np
from pyretis.orderparameter import OrderParameter  
from itertools import product

logger = logging.getLogger(__name__)  
logger.addHandler(logging.NullHandler()) 


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

    This class defines the distance between the center of mass of 2 groups of particles
    in the drug delivery example.

    Attributes
    ----------
    name : string
        A human readable name for the order parameter
    index : integer
        This selects the particle to use for the order parameter.
    periodic : boolean
        This determines if periodic boundaries should be applied to
        the position or not.
    """

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

        Parameters
        ----------
        name : string
            The name for the order parameter
        index : tuple of ints
            This is the indices of the atom we will use the position of.
        periodic : boolean, optional
            This determines if periodic boundary conditions should be
            applied to the position.
        """
        super().__init__(description='Contact map')

        self.idx1=np.array([1086, 1098, 1103, 1110, 1123, 1126, 1129], dtype=np.int16)-1
        self.idx2=np.array([118,151,182,215,246,279,499,532,563,596,627,660], dtype=np.int16)-1


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

        Here, the order parameter is just the distance between the center of mass of the two groups of atoms defined above. 

        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.
        """

        r0=0.4  ; d0=0.25 ; nn=2  ; mm=4
        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)
        orderp=[0]                                                                            
        sij=[]
        for cnt1,in1 in enumerate(self.idx1):
            ia1=[] ; ia1.append(in1)
            for dist in md.compute_distances(trj[:],product(ia1,self.idx2)):
                for rij in dist:
                    rij=rij-d0
                    if rij < 0:
                        sij.append(1)
                    else:
                        sij.append(((1-(rij/r0)**nn)/(1-(rij/r0)**mm)))
        orderp[0]=sum(sij)

            # remove file_gro:
        if os.path.isfile(file_gro):
             logger.debug('Removing: %s', file_gro)
             os.remove(file_gro)
        print(orderp)
        return orderp[0]