molshape.pl

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


# *** Print help
# ***
if ($#ARGV < 2) {
  print("\nMOLSHAPE <G03-Cube-File>\n         <Atom number for start slice>\n         <Atom number for end slice>\n         <Axis definition (xyz|xzy|zyx|zxy), default = xyz>\n         <Electron density cut-off, default = 0.002>\n\n");
  exit;
}

$axis = lc($ARGV[3]);
if ($axis ne 'xyz'  &&  $axis ne 'xzy'  &&
    $axis ne 'yxz'  &&  $axis ne 'yzx'  &&
    $axis ne 'zxy'  &&  $axis ne 'zyx'  &&
    $axis ne 'xy'   &&  $axis ne 'xz'   &&
    $axis ne 'yx'   &&  $axis ne 'yz'   &&
    $axis ne 'zx'   &&  $axis ne 'zy') {
  $axis = 'xyz';            # Standard axis system
}

$border = $ARGV[4];
if ($border eq '') {
  $border = 0.002;          # Standard cut-off value
}

$bohr2angstrom = 0.5291772; # Conversion factor Bohr to Angstrom
$aaxis = substr($axis,0,1); # a-axis
$baxis = substr($axis,1,1); # b-axis
$caxis = substr($axis,2,1); # c-axis
$threedim = 1;              # auxiliary variable whether MOLSHAPE is in 3D mode (=1) or 2D mode (=0)
if ($axis eq 'xy') { $caxis = 'z'; $threedim = 0; }
if ($axis eq 'xz') { $caxis = 'y'; $threedim = 0; }
if ($axis eq 'yx') { $caxis = 'z'; $threedim = 0; }
if ($axis eq 'yz') { $caxis = 'x'; $threedim = 0; }
if ($axis eq 'zx') { $caxis = 'y'; $threedim = 0; }
if ($axis eq 'zy') { $caxis = 'x'; $threedim = 0; }

# *** Print title
# ***
$VERSION = 0.4;
$TITLE = sprintf("MOLSHAPE (Version %01.1f)",$VERSION);
print ("\n"."-" x length($TITLE)."\n");
print ($TITLE."\n");
print ("-" x length($TITLE)."\n\n");


# *** Open file
# ***
open(INPUTFILE, $ARGV[0]) or die("\nERROR: Cannot open file '$ARGV[0]'\n\n");
@inputlines = <INPUTFILE>;
close(INPUTFILE);


# *** Read header
# ***
foreach (0..5) { push(@header, $inputlines[$_]); }


# *** Read number of atoms and coordinate origin
# ***
$header[2]=~/\ *(\d*)\ *(-?\d*.\d*)\ *(-?\d*.\d*)\ *(-?\d*.\d*)/;
$atoms = $1;
if (substr($axis,0,2) eq 'xy') { $aorigin = $2; $borigin = $3; $corigin = $4; }
if (substr($axis,0,2) eq 'xz') { $aorigin = $2; $borigin = $4; $corigin = $3; }
if (substr($axis,0,2) eq 'yx') { $aorigin = $3; $borigin = $2; $corigin = $4; }
if (substr($axis,0,2) eq 'yz') { $aorigin = $3; $borigin = $4; $corigin = $2; }
if (substr($axis,0,2) eq 'zx') { $aorigin = $4; $borigin = $2; $corigin = $3; }
if (substr($axis,0,2) eq 'zy') { $aorigin = $4; $borigin = $3; $corigin = $2; }


# *** Read number of voxels and voxel lengths
# ***
$header[3]=~/\ *(\d*)\ *(\d*.\d*)\ */;
  $aavoxels = $1; # auxiliary variable for cube field
  if ($aaxis eq 'x') { $avoxels = $1; $avoxellength = $2; }
  if ($baxis eq 'x') { $bvoxels = $1; $bvoxellength = $2; }
  if ($caxis eq 'x') { $cvoxels = $1; $cvoxellength = $2; }
$header[4]=~/\ *(\d*)\ *\d*.\d*\ *(\d*.\d*)\ */;
  $bbvoxels = $1; # auxiliary variable for cube field
  if ($aaxis eq 'y') { $avoxels = $1; $avoxellength = $2; }
  if ($baxis eq 'y') { $bvoxels = $1; $bvoxellength = $2; }
  if ($caxis eq 'y') { $cvoxels = $1; $cvoxellength = $2; }
$header[5]=~/\ *(\d*)\ *\d*.\d*\ *\d*.\d*\ *(\d*.\d*)/;
  $ccvoxels = $1; # auxiliary variable for cube field
  if ($aaxis eq 'z') { $avoxels = $1; $avoxellength = $2; }
  if ($baxis eq 'z') { $bvoxels = $1; $bvoxellength = $2; }
  if ($caxis eq 'z') { $cvoxels = $1; $cvoxellength = $2; }


# *** Print out the read data
# ***
if ($threedim) {
  print ("Axis definition         : $aaxis (main axis), $baxis (secondary axis), $caxis (secondary axis)\n");
}
else {
  print ("Axis definition         : $aaxis (main axis), $baxis (secondary axis)\n");
}
print ("Electron density cut-off: $border electrons/Bohr^3\n");
print ("Number of atoms         : $atoms\n\n");
print ("Number/length $aaxis-voxels  : $avoxels/$avoxellength Bohr\n");
print ("Number/length $baxis-voxels  : $bvoxels/$bvoxellength Bohr\n");
print ("Number/length $caxis-voxels  : $cvoxels/$cvoxellength Bohr\n");
$voxelarea = $bvoxellength*$cvoxellength*$bohr2angstrom**2;
printf ("Voxel area ($baxis,$caxis)        : %.5f Bohr^2  =  %.5f A^2\n\n",$voxelarea/($bohr2angstrom**2),$voxelarea);


# *** Start of the cube definition
# ***
$cubestart=$atoms+2+1+3; # *** $atoms atom coordinates + TWO comment lines + ONE origin definition + THREE voxel definitions


# *** Read cube definition
# ***
foreach $l ($cubestart..$#inputlines) # iterate through cube definition
{
  foreach (split(/\s+/,substr($inputlines[$l],2))) {  # TWO leading spaces removed, then split by whitespace
    push(@cube, $_);
  }
}
if ($#cube-($avoxels*$bvoxels*$cvoxels-1)!=0) {
  die("ERROR: The read cube definition does not match the cube dimensions!\n");
}
else {
  print ("OK: The read cube definition matches the cube dimensions!\n\n");
}

# *** Convert cube definition into cube array (a,b,c)
# ***
foreach $ia (1..$aavoxels) {
  foreach $ib (1..$bbvoxels) {
    foreach $ic (1..$ccvoxels) {
      # DENSITY = $cubearray[$xvox][$yvox][$zvox][0];
      # FLAG = $cubearray[$xvox][$yvox][$zvox][1];
      if (substr($axis,0,2) eq 'xy') { $cubearray[$ia][$ib][$ic]=[shift(@cube),0]; }
      if (substr($axis,0,2) eq 'xz') { $cubearray[$ia][$ic][$ib]=[shift(@cube),0]; }
      if (substr($axis,0,2) eq 'yx') { $cubearray[$ib][$ia][$ic]=[shift(@cube),0]; }
      if (substr($axis,0,2) eq 'yz') { $cubearray[$ib][$ic][$ia]=[shift(@cube),0]; }
      if (substr($axis,0,2) eq 'zx') { $cubearray[$ic][$ia][$ib]=[shift(@cube),0]; }
      if (substr($axis,0,2) eq 'zy') { $cubearray[$ic][$ib][$ia]=[shift(@cube),0]; }
    }
  }
}


# *** Create empty MolShape array (a,b,c)
# ***
foreach $ia (1..$aavoxels) {
  foreach $ib (1..$bbvoxels) {
    foreach $ic (1..$ccvoxels) {
      $molshapearray[$ia][$ib][$ic]=0;
    }
  }
}


# *** Assumed cavity/volume center (b,c)
# ***
$bvoxelstart = ($bvoxels-($bvoxels%2))/2;
$cvoxelstart = ($cvoxels-($cvoxels%2))/2;
print ("Cube-Zentrum ($baxis,$caxis)      : ($bvoxelstart,$cvoxelstart)\n");


# *** Determine a-slices
# ***
$inputlines[$ARGV[1]+2+1+3-1] =~ /\ *\d*\ *\d*.\d*\ *(-?\d*.\d*)\ *(-?\d*.\d*)\ *(-?\d*.\d*)/;
if ($aaxis eq 'x') { $astartcoord = $1; }
if ($aaxis eq 'y') { $astartcoord = $2; }
if ($aaxis eq 'z') { $astartcoord = $3; }
$astartcoordapprox = $aorigin;
$astartslice = 0;
while ($astartcoordapprox<=$astartcoord) { # search start coordinate
  $astartcoordapprox += $avoxellength;
  $astartslice++;
}

$inputlines[$ARGV[2]+2+1+3-1]=~/\ *\d*\ *\d*.\d*\ *(-?\d*.\d*)\ *(-?\d*.\d*)\ *(-?\d*.\d*)/;
if ($aaxis eq 'x') { $aendcoord = $1; }
if ($aaxis eq 'y') { $aendcoord = $2; }
if ($aaxis eq 'z') { $aendcoord = $3; }
$aendcoordapprox = $aorigin;
$aendslice = 0;
while ($aendcoordapprox<=$aendcoord) { # search end coordinate
  $aendcoordapprox += $avoxellength;
  $aendslice++;
}

if ($aendslice<$astartslice) { # swap start and end if end coordinates < start coordinates
  ($aendslice,$astartslice) = ($astartslice,$aendslice);
  ($aendcoord,$astartcoord) = ($astartcoord,$aendcoord);
  ($ARGV[1],$ARGV[2]) = ($ARGV[2],$ARGV[1])
}

print ("first $aaxis-voxel           : $astartslice (atom number $ARGV[1])\n");
print ("last $aaxis-voxel            : $aendslice (atom number $ARGV[2])\n\n");


# *** Volume or cavity mode?
# ***
$volume = undef; # auxiliary variable for volume (=1) or cavity mode (=0)
$avoxelmiddle = int(($aendslice+$astartslice)/2);
if ($cubearray[$avoxelmiddle][$bvoxelstart][$cvoxelstart][0] > $border) {
  $volume = 1; # volume mode
  print ("MolShape mode           : Volume mode\n\n");
}
else {
  $volume = 0; # cavity mode
  print ("MolShape mode           : Cavity mode\n\n");
}


# *** Create profile
# ***
$cavevolume = 0;
$output = $ARGV[0];
if ($volume) { # volume mode
  $outputinfo = "_volume_".$axis."_".$border;
}
else {         # cavity mode
  $outputinfo = "_cavity_".$axis."_".$border;
}
$output =~ s/\.cube|\.cub/$outputinfo.csv/;
open (OUTPUTFILE, ">$output") or die("ERROR: Cannot create file '$output'!\n\n");
print OUTPUTFILE ("slice, coord_rel, coord_abs, slicearea, integration, equiv_radius\n");
foreach ($astartslice..$aendslice) {
  $voxelcount=0;
  if ($volume) {
    BoundaryFillVolume($bvoxelstart,$cvoxelstart,$_); # volume mode
  }
  else {
    BoundaryFillCavity($bvoxelstart,$cvoxelstart,$_); # cavity mode
  }
  $caveslicearea = $voxelcount*$voxelarea;  # area of the current slice (in A^2)
  $cavevolume += $caveslicearea*$avoxellength*$bohr2angstrom; # molecule/cavity volume (in A^3)
  $equivradius = sqrt($caveslicearea/pi);   # equivalent radius of the current slice (in A)
  printf OUTPUTFILE ("%u, %.4f, %.4f, %.2f, %.2f, %.3f\n", $_, ($_-$astartslice)*$avoxellength*$bohr2angstrom, ($aorigin+$_*$avoxellength)*$bohr2angstrom, $caveslicearea, $cavevolume, $equivradius);
 }
close(OUTPUTFILE);
print ("OK: The file '$output' was created successfully.\n");


# *** Write new cube file
# ***
$output = $ARGV[0];
$output =~ s/\.cube|\.cub/$outputinfo.cube/;
open(OUTPUTFILE, ">$output") or die("ERROR: Cannot create file '$output'!\n\n");
foreach (0..$cubestart-1) {
  print OUTPUTFILE $inputlines[$_];
}
foreach $ia (1..$aavoxels) {
  foreach $ib (1..$bbvoxels) {
    foreach $ic (1..$ccvoxels) {
      if (substr($axis,0,2) eq 'xy') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ia][$ib][$ic]); }
      if (substr($axis,0,2) eq 'xz') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ia][$ic][$ib]); }
      if (substr($axis,0,2) eq 'yx') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ib][$ia][$ic]); }
      if (substr($axis,0,2) eq 'yz') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ib][$ic][$ia]); }
      if (substr($axis,0,2) eq 'zx') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ic][$ia][$ib]); }
      if (substr($axis,0,2) eq 'zy') { printf OUTPUTFILE ("  %.5E",$molshapearray[$ic][$ib][$ia]); }
      if ($ic%6 == 0) { print OUTPUTFILE ("\n"); }
    }
    print OUTPUTFILE ("\n");
  }
}
close(OUTPUTFILE);
print ("OK: The file '$output' was created successfully.\n\n");




# *** Define BoundaryFill functions
# ***
sub BoundaryFillVolume { # (y,z,x)
  local $bvox = shift(@_);
  local $cvox = shift(@_);
  local $avox = shift(@_);
  if ($cubearray[$avox][$bvox][$cvox][0]>$border and $cubearray[$avox][$bvox][$cvox][1]==0 and $bvox>=1 and $bvox<=$bvoxels and $cvox>=1 and $cvox<=$cvoxels) {
    $cubearray[$avox][$bvox][$cvox][1]=1;  # voxel is flagged
    $voxelcount++; # "volume voxels" + 1
    $molshapearray[$avox][$bvox][$cvox]=1; # electron density of the voxel > $border
    BoundaryFillVolume($bvox+1,$cvox,$avox);
    BoundaryFillVolume($bvox-1,$cvox,$avox);
    if ($threedim) { BoundaryFillVolume($bvox,$cvox+1,$avox) };
    if ($threedim) { BoundaryFillVolume($bvox,$cvox-1,$avox) };
  }
 return;
}

sub BoundaryFillCavity { # (y,z,x)
  local $bvox = shift(@_);
  local $cvox = shift(@_);
  local $avox = shift(@_);
  if ($cubearray[$avox][$bvox][$cvox][0]<$border and $cubearray[$avox][$bvox][$cvox][1]==0 and $bvox>=1 and $bvox<=$bvoxels and $cvox>=1 and $cvox<=$cvoxels) {
    $cubearray[$avox][$bvox][$cvox][1]=1;  # voxel is flagged
    $voxelcount++; # "cavity voxels" + 1
    $molshapearray[$avox][$bvox][$cvox]=1; # electron density of the voxel < $border
    BoundaryFillCavity($bvox+1,$cvox,$avox);
    BoundaryFillCavity($bvox-1,$cvox,$avox);
    if ($threedim) { BoundaryFillCavity($bvox,$cvox+1,$avox) };
    if ($threedim) { BoundaryFillCavity($bvox,$cvox-1,$avox) };
  }
 return;
}