#!/bin/bash

# Author: David Small, MGCF, May 2023

# usage: charge_centers_from_mae file.mae
# file.mae must contain (force-field) partial charges
#	one way to get this is to use desmond system builder, then export that entry as mae
# the script uses the r_m_charge1 column from file.mae, because that's what compute_dipole_charged.py uses


m_atom_line=$(grep -n 'm_atom\[' $1 | awk '{ split($1,arr,":"); print arr[1]; }')

tail -n +$m_atom_line $1 | awk 'BEGIN { 

 readData = 0;
 n_charge = 0.0; p_charge = 0.0;
 for (i=1; i<=6; i++) nr[i] = 0.0;

}

{

 if (readData)
 {
  if ($1 == ":::") exit;
  else
  {
	n=split($0,arr,"\"");

	count = 1;

	for (i = 1; i <= n; i++)
	{
        	if ((i % 2) == 0)
        	{
                	lineCols[count] = "";
 	        	count++;
        	}
		else
		{
                	nn = split(arr[i], arrarr, " ");
                	for (ii = 1; ii <= nn; ii++)
                	{
                        	lineCols[count] = arrarr[ii];
 	                	count++;
                	}
        	}
	}

	# print count;


    chg = lineCols[qcol];
    if (chg > 0.0)
    {
	p_charge += chg;
	j = 3;
    }
    else
    {
	n_charge += chg;
        j = 0;
    }
    for (i=1; i<4; i++) nr[j+i] += chg * lineCols[rcol[i]];
  }
 }
 else
 {
  if ($1 == "r_m_x_coord") rcol[1] = NR-1;
  if ($1 == "r_m_y_coord") rcol[2] = NR-1;
  if ($1 == "r_m_z_coord") rcol[3] = NR-1;
  if ($1 == "r_m_charge1") qcol = NR-1;

  if ($1 == ":::") readData = 1;
 }

}

END {

 for (i=1; i<=3; i++)
 {
  if (n_charge < 0.0) nr[i] /= n_charge;
  if (p_charge > 0.0) nr[i+3] /= p_charge;
 }

 print "Total negative charge:",n_charge;
 print "Center of negative charge:",nr[1],nr[2],nr[3];
 print "Total positive charge:",p_charge;
 print "Center of positive charge:",nr[4],nr[5],nr[6];

}'