Hello Monks! Please save me from my impending insanity! I am desperate to finish this assignment and just can't get it to work properly despite spending hours on it. The assignment involves writing a progras that uses a Metropolis Monte Carlo algorithm to generate random configurations of noble gas atoms in a hard box... overall it should result a lowering of energy by moving the configurations of the atoms in space. I have the feeling that the problem lies with my energy calculations of both the initial/previous energy ($prevE) and the new energy ($newE). For some reason, my previous and new energies are often identical despite moving am atom in space (confirmed this my printing $prevE and $newE)... I have no idea why this is happening as obviously moving atoms should have a direct impact on energy... Any suggestions will be greatly appreciated. PS. sorry for my poor quality coding, but I am a chemist attempting to learn perl as part of my honours course. 

#!/usr/bin/perl -w

$loops = 5000;             # no. of mc loops
@box = (10, 10, 10);     # box dimensions
$k = 0.001987204118;     #kcal/mole/K
$T = 298;                 #K
$prevE = 0;
$atom_name = 'Ar';
$residue_name = 'NOB B';

open (OUTFILE, ">pdb_output.pdb") || die "could not open file";

# store initial coordinates (from pdb file)
# calculate initial energy of the system

{
    $atom_number = 0;

        while (<>) {  

# Find x, y, z coordinates and store in separate arrays (all x, all y,
+ all z)

            if ($_ =~ /^ATOM/) {
                @line = $_ =~ m/^(.....).(.....).(....).(...)..(....).
+...(........)(........)(........)/;
        
                $x = $line[5];
                $arrayx[$atom_number] = $x;
                         
                $y = $line[6];
                $arrayy[$atom_number] = $y;
        
                $z = $line[7];
                $arrayz[$atom_number] = $z;
        
                ++$atom_number;        
            }    
        }
}
    
# Monte Carlo algorithm

foreach (1..$loops) {

# Calculate initial/previous energy
    
    
    $prevE = total_energy ();
    
    print"-----------------\nPrevious energy = $prevE\n";
    
# Pick a random atom

    $atom = int(rand($atom_number));

    print "Atom number = $atom\n";
    
# Save inital atom coordinates into an array

    @prev_coord = ($arrayx[$atom], $arrayy[$atom], $arrayz[$atom]);

    print "Initial coordinates = @prev_coord\n";

# Rondomly move atom along each x, y and z axis

    $dx = rand (2) - 1;
    $dy = rand (2) - 1;
    $dz = rand (2) - 1;

    $arrayx[$atom] = $arrayx[$atom] + $dx;
    $arrayy[$atom] = $arrayy[$atom] + $dy;
    $arrayz[$atom] = $arrayz[$atom] + $dz;

# Make sure atom remains inside the box
     
     $arrayx[$atom] = in_box($arrayx[$atom], $box[0]);
    $arrayy[$atom] = in_box($arrayy[$atom], $box[1]);
    $arrayz[$atom] = in_box($arrayz[$atom], $box[2]);
    
    print "New coordinates = $arrayx[$atom], $arrayy[$atom], $arrayz[$
+atom]\n";
    
# Calculation of new energy

    $newE = total_energy ();
    
    print "New energy = $newE\n";
    
# Accept or reject algortithm

{
    if ($newE <= $prevE) {
    
    $decision = 'accept';
            
    } else {
    
        $R = rand 1.0;
        $B = exp(-($newE - $prevE) / ($k * $T));
    
        if ($R >= $B) {
    
        $decision = 'reject';
                    
        } else {
    
        $decision = 'accept';
        
        }
        
        print "R = $R >= B = $B;\n";
    }
    
    if ($decision eq 'accept') {
        
        $lowestE = $newE;
        
        for ($i = 0; $i <= $#arrayx; ++$i) {
        
            printf OUTFILE ("ATOM %5d %-4s %3s %4d %8.3f%8.3f%8.3f\n",
            $i+1, $atom_name, $residue_name, $i+1, $arrayx[$i], $array
+y[$i], $arrayz[$i]);
        
        }
        
        print OUTFILE "END\n";
        
        print "Decision = Accept\n";
        
        } else {
        
        $lowestE = $prevE;
        
        #restore previous coordinates

        $arrayx[$atom] = $prev_coord[0];
        $arrayy[$atom] = $prev_coord[1];
        $arrayz[$atom] = $prev_coord[2];
            
        print "Decision = Reject\n";
        }


        $decision = 0;
        print "Energy = $lowestE\n";

        
    }
    
    }    
    
close OUTFILE;        
    
    
#################    
## SUBROUTINES ##    
#################

# LJ energy calculation subroutine

sub vdw {

    my $r = $_[0];

    my $s = 2.780;
    my $e = 0.069;
       

       $LJ_E = 4*$e*(
                       ($s/$r)**12 - 
                       ($s/$r)**6
                       );
                       
    return $LJ_E;
}

# Keeping atom in box subroutine

sub in_box {

    my $coord = $_[0];
    my $box = $_[1];
    
    if ($coord > $box) {
        $coord = $coord - $box;
    } else {
        $coord = $coord;
    }
    
    if ($coord < 0) {
        $coord = $coord + $box;
    } else {
        $coord = $coord;
    }
        return $coord;    
}    

sub total_energy {
    
    $totalE = 0;

    foreach my $i (0..$#arrayx) {

        foreach my $j ($i + 1..$#arrayx) {
    
            $r = sqrt(
                    ($arrayx[$i] - $arrayx[$j])**2 + 
                    ($arrayy[$i] - $arrayy[$j])**2 + 
                    ($arrayz[$i] - $arrayz[$j])**2
                    );        
    
            $LJ_E = vdw($r);

            $totalE = $totalE + $LJ_E;
            
            return $totalE;
        }
    }
}

[download]
The pdb file that it reads from is: 
HEADER    noble
COMPND    noble gas
ATOM      1  Ar  NOB B   1      -2.252   0.390   2.441  1.00 11.60    
+       N
ATOM      2  Ar  NOB B   2       0.482  -0.812   2.262  1.00 10.61    
+       C
ATOM      3  Ar  NOB B   3       1.490  -1.001  -2.448  1.00 12.28    
+       N
ATOM      4  Ar  NOB B   4       2.430  -0.614   3.444  1.00 12.22    
+       C
ATOM      5  Ar  NOB B   5       3.948   0.553   1.512  1.00 11.87    
+       C
ATOM      6  Ar  NOB B   6      -3.947   0.986  -1.075  1.00 15.17    
+       N
ATOM      7  Ar  NOB B   7       0.344   1.001  -2.335  1.00  9.75    
+       C
ATOM      8  Ar  NOB B   8      -1.696   0.570  -3.444  1.00 11.08    
+       C
END

[download]
Thanks for your help, its so very greatly appreciated :)

In reply to Metropolis Monte Carlo by stellaparallax

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post, it's "PerlMonks-approved HTML":


  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.