#!/usr/bin/perl use Math::Trig; open(DATA,">data.brown") || die "data.brown could not be opened"; printf DATA ("# heated Brownian particle\n\n"); $cutoff = 1.36; $cutoff2 = $cutoff*$cutoff; # box $rhos = 1.59; $as = (4.0/$rhos)**(1.0/3.0); $R = 4.4; $ns = int($R/$as)+1; $rhol = 0.76; $al = (4.0/$rhol)**(1.0/3.0); $nl = 9; $xlo = -$nl*$al; $ylo = -$nl*$al; $zlo = -$nl*$al; $xhi = $nl*$al; $yhi = $nl*$al; $zhi = $nl*$al; # solid $atom = 0; for ($ix=-$ns;$ix<$ns;$ix++) { for ($iy=-$ns;$iy<$ns;$iy++) { for ($iz=-$ns;$iz<$ns;$iz++) { $x = $ix*$as; $y = $iy*$as; $z = $iz*$as; if ($x*$x+$y*$y+$z*$z<$R*$R) { $atom++; $at_type[$atom] = 2; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = ($ix+0.5)*$as; $y = ($iy+0.5)*$as; $z = $iz*$as; if ($x*$x+$y*$y+$z*$z<$R*$R) { $atom++; $at_type[$atom] = 2; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = $ix*$as; $y = ($iy+0.5)*$as; $z = ($iz+0.5)*$as; if ($x*$x+$y*$y+$z*$z<$R*$R) { $atom++; $at_type[$atom] = 2; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = ($ix+0.5)*$as; $y = $iy*$as; $z = ($iz+0.5)*$as; if ($x*$x+$y*$y+$z*$z<$R*$R) { $atom++; $at_type[$atom] = 2; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } } } } $N_s = $atom; $bond = 0; for ($i=1;$i<=$atom;$i++) { for ($j=$i+1;$j<=$atom;$j++) { $dr2 = ($at_x[$i]-$at_x[$j])*($at_x[$i]-$at_x[$j])+($at_y[$i]-$at_y[$j])*($at_y[$i]-$at_y[$j])+($at_z[$i]-$at_z[$j])*($at_z[$i]-$at_z[$j]); if ($dr2 < $cutoff2) { $bond++; $bond_at1[$bond] = $i; $bond_at2[$bond] = $j; } } } # liquid for ($ix=-$nl;$ix<$nl;$ix++) { for ($iy=-$nl;$iy<$nl;$iy++) { for ($iz=-$nl;$iz<$nl;$iz++) { $x = $ix*$al; $y = $iy*$al; $z = $iz*$al; if ($x*$x+$y*$y+$z*$z>($R+1.0)*($R+1.0)) { $atom++; $at_type[$atom] = 1; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = ($ix+0.5)*$al; $y = ($iy+0.5)*$al; $z = $iz*$al; if ($x*$x+$y*$y+$z*$z>($R+1.0)*($R+1.0)) { $atom++; $at_type[$atom] = 1; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = $ix*$al; $y = ($iy+0.5)*$al; $z = ($iz+0.5)*$al; if ($x*$x+$y*$y+$z*$z>($R+1.0)*($R+1.0)) { $atom++; $at_type[$atom] = 1; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } $x = ($ix+0.5)*$al; $y = $iy*$al; $z = ($iz+0.5)*$al; if ($x*$x+$y*$y+$z*$z>($R+1.0)*($R+1.0)) { $atom++; $at_type[$atom] = 1; $at_x[$atom] = $x; $at_y[$atom] = $y; $at_z[$atom] = $z; } } } } $N_at = $atom; $N_bond = $bond; print "N_at = $N_at N_s = $N_s N_bond = $N_bond\n"; printf DATA ("$N_at\t atoms\n"); printf DATA ("$N_bond\t bonds\n\n"); printf DATA ("2\t atom types\n"); printf DATA ("1\t bond 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 ("\nAtoms\n\n"); for ($atom=1;$atom<=$N_at;$atom++) { printf DATA ("$atom 0 $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 ("\n"); close (DATA);