#!/usr/bin/perl
use Math::Trig;

open(DATA,">data.cnt") || die
  "data.cnt could not be opened";
open(DUMP,">dump0.cnt.xyz") || die
  "dump0.cnt.xyz could not be opened";
printf DATA ("# capillary filling of a SWCNT with water\n\n");

# == parameters ==

# Masses
$mH = 1.008;
$mO = 16.00;
$mC = 12.01;

# Interactions
$epsH = 0.0;
$epsO = 0.1521;
$epsC = 0.086; 
$epsC1 = 0.00357;

$sigmaH = 0.0;
$sigmaO = 3.1507;
$sigmaC = 3.400;

# carbon
$typeC = 6;
$qC = 0.0;

# TIP3P water
$dOH = 0.9572;
$HOH = 104.52*pi/180.;
$qO = -0.834;
$qH = 0.417;
$typeO = 8;
$typeH = 1;

# graphene and CNT
$d = 1.42;  # C-C distance
$a1 = sqrt(3.0)*$d;
$a2 = 3.0*$d;

$na1 = 40;  # number of unit cells along the axis of the CNT
$na2 = 30;  # (na2,na2) armchair nanotube

$radius = 3.0*$na2*$d/(2.0*pi);  # radius in angstroms
$diameter = 2.0*$radius;  # diameter in angstroms
$r_eff = $radius-0.25*($sigmaO + $sigmaC);
$length = $na1*$a1;  # length in angstroms
print "radius = $radius A; r_eff = $r_eff A; length = $length A\n";

# tank dimensions
$rho = 0.0325;
$Nmol_CNT = $rho*(pi*($r_eff)**(2.0))*$length;
print "Nmol_CNT ~ $Nmol_CNT\n";

$nax = int(50.0/$a1)+1;
$nay = int(50.0/$a2)+1;

$h = 100.0;

# water initial sc lattice
$a = $rho**(-1.0/3.0);
$nz = int(($h-2.0*sqrt($sigmaO * $sigmaC)+$sigmaO)/$a);
$z0 = -$nz*$a;

$xlo = -$nax*$a1;
$ylo = -$nay*$a2;
$zlo = -$h-22.0;
$xhi = $nax*$a1;
$yhi = $nay*$a2;
$zhi = $length+7.0;

# water initial sc lattice
$nx = int(($xhi-$xlo)/$a);
$ny = int(($yhi-$ylo)/$a);
print "water initial sc lattice: nx = $nx  ny = $ny  nz = $nz\n";


# == carbon nanotube ==

$atom = 0;
$molecule = 0;
for ($ia1=0;$ia1<$na1;$ia1++) {
  for ($ia2=0;$ia2<$na2;$ia2++) {
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC;
    $at_q[$atom] = $qC;
    $at_z[$atom] = $ia1*$a1;
    $theta = $ia2*$a2/$radius;
    $at_x[$atom] = $radius*cos($theta);
    $at_y[$atom] = $radius*sin($theta);
    
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC;
    $at_q[$atom] = $qC;
    $at_z[$atom] = $ia1*$a1;
    $theta = ($ia2*$a2 + 2.0*$d)/$radius;
    $at_x[$atom] = $radius*cos($theta);
    $at_y[$atom] = $radius*sin($theta);
    
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC;
    $at_q[$atom] = $qC;
    $at_z[$atom] = ($ia1+0.5)*$a1;
    $theta = ($ia2*$a2 + 0.5*$d)/$radius;
    $at_x[$atom] = $radius*cos($theta);
    $at_y[$atom] = $radius*sin($theta);
    
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC;
    $at_q[$atom] = $qC;
    $at_z[$atom] = ($ia1+0.5)*$a1;
    $theta = ($ia2*$a2 + 1.5*$d)/$radius;
    $at_x[$atom] = $radius*cos($theta);
    $at_y[$atom] = $radius*sin($theta);
  }
}


# == graphene sheets ==

$dist2test = ($radius + (0.99)*$d)*($radius + (0.99)*$d);
$dist2testplug = ($radius - (0.99)*$d)*($radius - (0.99)*$d);
$z = 0.0;
for ($ia1=-$nax;$ia1<$nax;$ia1++) {
  for ($ia2=-$nay;$ia2<$nay;$ia2++) {
    $x = $ia1*$a1;
    $y = $ia2*$a2;
    $dist2 = $x*$x+$y*$y;
    if (($dist2 > $dist2test)||($dist2 < $dist2testplug)) {
      $atom++;
      $at_molecule[$atom] = $molecule;
      #$at_type[$atom] = $typeC;
      $at_type[$atom] = ($dist2>$dist2test)?$typeC:$typeC+1;
      $at_q[$atom] = $qC;
      $at_x[$atom] = $x;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
    }
    
    #$x = $ia1*$a1;
    $y = $ia2*$a2 + 2.0*$d;
    $dist2 = $x*$x+$y*$y;
    if (($dist2 > $dist2test)||($dist2 < $dist2testplug)) {
      $atom++;
      $at_molecule[$atom] = $molecule;
      #$at_type[$atom] = $typeC;
      $at_type[$atom] = ($dist2>$dist2test)?$typeC:$typeC+1;
      $at_q[$atom] = $qC;
      $at_x[$atom] = $x;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
    }
    
    $x = ($ia1+0.5)*$a1;
    $y = $ia2*$a2 + 0.5*$d;
    $dist2 = $x*$x+$y*$y;
    if (($dist2 > $dist2test)||($dist2 < $dist2testplug)) {
      $atom++;
      $at_molecule[$atom] = $molecule;
      #$at_type[$atom] = $typeC;
      $at_type[$atom] = ($dist2>$dist2test)?$typeC:$typeC+1;
      $at_q[$atom] = $qC;
      $at_x[$atom] = $x;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
    }
    
    #$x = $ia1*$a1 + 0.5*sqrt(3.0)*$d;
    $y = $ia2*$a2 + 1.5*$d;
    $dist2 = $x*$x+$y*$y;
    if (($dist2 > $dist2test)||($dist2 < $dist2testplug)) {
      $atom++;
      $at_molecule[$atom] = $molecule;
      #$at_type[$atom] = $typeC;
      $at_type[$atom] = ($dist2>$dist2test)?$typeC:$typeC+1;
      $at_q[$atom] = $qC;
      $at_x[$atom] = $x;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
    }
  }
}

$N_piston = 0;
$z = -$h-15.0;
for ($ia1=-$nax;$ia1<$nax;$ia1++) {
  for ($ia2=-$nay;$ia2<$nay;$ia2++) {
    $x = $ia1*$a1;
    $y = $ia2*$a2;
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC-1;
    $at_q[$atom] = $qC;
    $at_x[$atom] = $x;
    $at_y[$atom] = $y;
    $at_z[$atom] = $z;
    $N_piston ++;
    
    #$x = $ia1*$a1;
    $y = $ia2*$a2 + 2.0*$d;
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC-1;
    $at_q[$atom] = $qC;
    $at_x[$atom] = $x;
    $at_y[$atom] = $y;
    $at_z[$atom] = $z;
    $N_piston ++;
    
    $x = ($ia1+0.5)*$a1;
    $y = $ia2*$a2 + 0.5*$d;
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC-1;
    $at_q[$atom] = $qC;
    $at_x[$atom] = $x;
    $at_y[$atom] = $y;
    $at_z[$atom] = $z;
    $N_piston ++;
    
    #$x = $ia1*$a1 + 0.5*sqrt(3.0)*$d;
    $y = $ia2*$a2 + 1.5*$d;
    $atom++;
    $at_molecule[$atom] = $molecule;
    $at_type[$atom] = $typeC-1;
    $at_q[$atom] = $qC;
    $at_x[$atom] = $x;
    $at_y[$atom] = $y;
    $at_z[$atom] = $z;
    $N_piston ++;
    
  }
}


# == water ==

$molecule = 0;
$bond = 0;
$angle = 0;

for ($iz=0;$iz<$nz;$iz++) {
  $z = $z0 + $iz*$a;
  for ($iy=0;$iy<$ny;$iy++) {
    $y = $ylo + $iy*$a;
    for ($ix=0;$ix<$nx;$ix++) {
      $x = $xlo + $ix*$a;
      $atom++;
      $molecule++;
      $at_molecule[$atom] = $molecule;
      $at_type[$atom] = $typeO;
      $at_q[$atom] = $qO;
      $at_x[$atom] = $x;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
      
      $atom++;
      $at_molecule[$atom] = $molecule;
      $at_type[$atom] = $typeH;
      $at_q[$atom] = $qH;
      $at_x[$atom] = $x+$dOH;
      $at_y[$atom] = $y;
      $at_z[$atom] = $z;
      
      $atom++;
      $at_molecule[$atom] = $molecule;
      $at_type[$atom] = $typeH;
      $at_q[$atom] = $qH;
      $at_x[$atom] = $x+$dOH*cos($HOH);
      $at_y[$atom] = $y+$dOH*sin($HOH);
      $at_z[$atom] = $z;
      
      $bond++;
      $bond_at1[$bond] = $atom - 2;
      $bond_at2[$bond] = $atom - 1;
      $bond++;
      $bond_at1[$bond] = $atom - 2;
      $bond_at2[$bond] = $atom;
      $angle++;
      $angle_at0[$angle] = $atom - 2;
      $angle_at1[$angle] = $atom - 1;
      $angle_at2[$angle] = $atom;
    }
  }
}

$N_at = $atom;
$N_bond = $bond;
$N_angle = $angle;
print "N_at = $N_at  N_mol = $molecule\n";

printf DUMP ("$N_at\nAtoms\n");

printf DATA ("$N_at\t atoms\n");
printf DATA ("$N_bond\t bonds\n");
printf DATA ("$N_angle\t angles\n\n");

printf DATA ("8\t atom types\n");
printf DATA ("1\t bond types\n");
printf DATA ("1\t angle types\n\n");

printf DATA ("$xlo $xhi xlo xhi\n");
printf DATA ("$ylo $yhi ylo yhi\n");
printf DATA ("$zlo $zhi zlo zhi\n\n");

printf DATA ("Masses\n\n");
for ($k=1;$k<=8;$k++) {
  $m = $mC;
  $m = $mH if ($k==$typeH);
  $m = $mO if ($k==$typeO);
  printf DATA ("$k\t $m\n");
}

printf DATA ("\nAtoms\n\n");

for ($atom=1;$atom<=$N_at;$atom++) {
  printf DATA ("$atom $at_molecule[$atom] $at_type[$atom] $at_q[$atom] $at_x[$atom] $at_y[$atom] $at_z[$atom]\n");
  printf DUMP ("$at_type[$atom] $at_x[$atom] $at_y[$atom] $at_z[$atom]\n");
}

printf DATA ("\nBonds\n\n");

for ($bond=1;$bond<=$N_bond;$bond++) {
  printf DATA ("$bond 1 $bond_at1[$bond] $bond_at2[$bond]\n");
}

printf DATA ("\nAngles\n\n");

for ($angle=1;$angle<=$N_angle;$angle++) {
  printf DATA ("$angle 1 $angle_at1[$angle] $angle_at0[$angle] $angle_at2[$angle]\n");
}

printf DATA ("\n");

close (DATA);
close (DUMP);