//methods for Verlet class. The constructors set the data members of the Verlet class. The initialize() function prepares the system to be time stepped by the integrate() function which contains the main molecular dynamics step loop.
#include "Verlet/Verlet.h"
#include "Verlet/Force.h"
#include <iostream.h>
#include <stdlib.h>
#include <math.h>


Verlet::Verlet() {
    cout << "enter the initial temperature.\n initial temperature = ";
    cin >> m_temperature;
    cout << "enter the number of particles.\n number of particles = ";
    cin >> m_numberOfParticles;
    cout << "enter the density.\n density = ";
    cin >> m_density;
    cout << "enter the time step.\n time step = ";
    cin >> m_timeStep;
    cout << "enter the number of time steps.\n number of time steps = ";
    cin >> m_numberOfTimeSteps;
}

Verlet::Verlet( const double temperature, const int numberOfParticles, const double density, const double timeStep, const int numberOfTimeSteps ) {
    m_temperature = temperature;
    m_numberOfParticles = numberOfParticles;
    m_density = density;
    m_timeStep = timeStep;
    m_numberOfTimeSteps = numberOfTimeSteps;
}

Verlet::~Verlet () {}

void Verlet::Initialize () {
  
  //edge length of entire volume of particles.
  m_bigEdge = pow ( ((double)(m_numberOfParticles))/m_density, 1./3. );

  //edge length of the little cube.
  m_littleEdge = pow ( m_density, -1./3. );

  //cube root of the number of particles in the volume.
  int cubeRootN = (int)(rint(pow ( ((double)(m_numberOfParticles)), 1./3. )));

  //center of mass velocities.
  double cmVx = 0.0;
  double cmVy = 0.0;
  double cmVz = 0.0;
  
  //sum of velocities squared
  double sumVsquared = 0.0; 

  for ( int zIter = 0; zIter<cubeRootN; zIter++ ) {
     for ( int yIter = 0; yIter<cubeRootN; yIter++ ) {
        for ( int xIter = 0; xIter<cubeRootN; xIter++ ) {

	   // initialize particle positions and velocities.
	   Particle aParticle( -m_bigEdge/2. + m_littleEdge/2. + ((double)(xIter))*m_littleEdge, -m_bigEdge/2. + m_littleEdge/2. + ((double)(yIter))*m_littleEdge, -m_bigEdge/2. + m_littleEdge/2. + ((double)(zIter))*m_littleEdge, ((double)(rand()))/((double)(RAND_MAX))-0.5, ((double)(rand()))/((double)(RAND_MAX))-0.5, ((double)(rand()))/((double)(RAND_MAX))-0.5 );
	   m_particleVector.push_back( aParticle );

	   cmVx += aParticle.Getvx();
	   cmVy += aParticle.Getvy();
	   cmVz += aParticle.Getvz();

	}
     }
  }

  cmVx /= m_numberOfParticles;
  cmVy /= m_numberOfParticles;
  cmVz /= m_numberOfParticles;

  // for debugging purposes.
  //cout << "big edge = " << m_bigEdge << endl << "little edge = " << m_littleEdge << endl << "cubeRootN = " << cubeRootN << endl << "CM vx = " << cmVx << endl << "CM vy = " << cmVy << endl << "CM vz = " << cmVz << endl;

  for ( vector<Particle>::iterator iter = m_particleVector.begin(); iter != m_particleVector.end(); iter++ ) {
     
    iter->Setvx( (iter->Getvx()-cmVx) );
    iter->Setvy( (iter->Getvy()-cmVy) );
    iter->Setvz( (iter->Getvz()-cmVz) );

    sumVsquared += iter->Getvx()*iter->Getvx() + iter->Getvy()*iter->Getvy() + iter->Getvz()*iter->Getvz();
  
  }

  //scale factor for velocities to ensure that total kinetic energy is zero.
  double scaleFactor = sqrt ( 3*((double)(m_numberOfParticles))*m_temperature/sumVsquared );

  for ( vector<Particle>::iterator iter = m_particleVector.begin(); iter != m_particleVector.end(); iter++ ) {
    
    iter->Setvx( iter->Getvx()*scaleFactor );
    iter->Setvy( iter->Getvy()*scaleFactor );
    iter->Setvz( iter->Getvz()*scaleFactor );

    iter->SetPrevx( iter->Getx() - iter->Getvx()*m_timeStep );
    iter->SetPrevy( iter->Gety() - iter->Getvy()*m_timeStep );
    iter->SetPrevz( iter->Getz() - iter->Getvz()*m_timeStep );

    //for debugging the particles initial conditions.
    //cout << "x = " << iter->Getx() << ", y = " << iter->Gety() << ", z = " << iter->Getz() << ", vx = " << iter->Getvx() << ", vy = " << iter->Getvy() << ", vz = " << iter->Getvz() << endl;

  } 
  
  //This block is used to debug the conservation of the intialization.
  /*
  cmVx = 0.0;
  cmVy = 0.0;
  cmVz = 0.0;
  sumVsquared = 0.0;
  
  for ( vector<Particle>::iterator iter = m_particleVector.begin(); iter != m_particleVector.end(); iter++ ) {
    
    cmVx += iter->Getvx();
    cmVy += iter->Getvy();
    cmVz += iter->Getvz();

    sumVsquared += (iter->Getvx())*(iter->Getvx()) + (iter->Getvy())*(iter->Getvy()) + (iter->Getvz())*(iter->Getvz());

  }

  cout << "cm vx = " << cmVx << endl << "cm vy = " << cmVy << endl << "cm vz = " << cmVz << endl << "sum v squared = " << sumVsquared << endl;
  */

}


void Verlet::Integrate ( vector<functionPtr> externalFunctions ) {
   
   //instantiate force function object.
   Force theForce;
 
   //local variables.
   double separationx = 0.0;
   double separationy = 0.0;
   double separationz = 0.0;
   double separation = 0.0;
   double xForce = 0.0;
   double yForce = 0.0;
   double zForce = 0.0;

   //time stepping loop.
   for ( int timeStep = 0; timeStep<m_numberOfTimeSteps; timeStep++ ) {

      //set forces to zero for every particle at each time step.
      for ( vector<Particle>::iterator particle = m_particleVector.begin(); particle != m_particleVector.end(); particle++ ) {
	 particle->zeroForce();
      } 

      //force calculation loop.
      for ( vector<Particle>::iterator i = m_particleVector.begin(); i != m_particleVector.end()-1; i++ ) {
         for ( vector<Particle>::iterator j = i+1; j != m_particleVector.end(); j++ ) {
	  
	   //calculate interparticle spacings using periodic boundary condition.
	   separationx = (i->Getx()-j->Getx());
	   separationx -= m_bigEdge*rint(separationx/m_bigEdge);
  	   separationy = (i->Gety()-j->Gety());
	   separationy -= m_bigEdge*rint(separationy/m_bigEdge);
	   separationz = (i->Getz()-j->Getz());
	   separationz -= m_bigEdge*rint(separationz/m_bigEdge);

	   separation = sqrt( (separationx*separationx) + (separationy*separationy) + (separationz*separationz) );


	   //for debugging the separation calculation.
	   /*
	     cout << "separation x = " << separationx << ", separation y = " << separationy << ", separation z = " << separationz << ", full sep = " << separation << endl;
	*/
	
	   xForce = theForce( separationx, separation );
	   yForce = theForce( separationy, separation );
	   zForce = theForce( separationz, separation );

	   i->updateXforce( xForce );
	   j->updateXforce( -xForce );
	   i->updateYforce( yForce );
	   j->updateYforce( -yForce );
	   i->updateZforce( zForce );
	   j->updateZforce( -zForce );

	 }
      }
   
      // used to debug force calculation
      /*
      for ( vector<Particle>::iterator particle = m_particleVector.begin(); particle != m_particleVector.end(); particle++ ) {
        cout << "x force = " << particle->GetXforce() << ", y force = " << particle->GetYforce() << ", z force = " << particle->GetZforce() << endl;
      }
      */

      //each particle integrates its own position and velocity individually.
      for ( vector<Particle>::iterator particle = m_particleVector.begin(); particle != m_particleVector.end(); particle++ ) {
	 particle->integrate ( m_timeStep );
      }

      //execute external functions.
      for ( vector<functionPtr>::iterator functionIter = externalFunctions.begin(); functionIter != externalFunctions.end(); functionIter++ ) {
	 functionIter->m_ptr2function( &m_temperature, &m_numberOfParticles, &m_density, &m_timeStep, &m_particleVector, &m_bigEdge, &m_littleEdge, &m_numberOfTimeSteps, &timeStep );
      }

   }

}
<!--'"</title></head>-->

<script type="text/javascript">
//OwnerIQ
var __oiq_pct = 50;
if( __oiq_pct>=100 || Math.floor(Math.random()*100/(100-__oiq_pct)) > 0 ) {
var _oiqq = _oiqq || [];
_oiqq.push(['oiq_addPageBrand','Lycos']);
_oiqq.push(['oiq_addPageCat','Internet > Websites']);
_oiqq.push(['oiq_addPageLifecycle','Intend']);
_oiqq.push(['oiq_doTag']);
(function() {
var oiq = document.createElement('script'); oiq.type = 'text/javascript'; oiq.async = true;
oiq.src = document.location.protocol + '//px.owneriq.net/stas/s/lycosn.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(oiq, s);
})();
}
//Google Analytics
var _gaq = _gaq || [];
_gaq.push(['_setAccount','UA-21402695-19']);
_gaq.push(['_setDomainName','tripod.com']);
_gaq.push(['_setCustomVar',1,'member_name','agentq314',3]);
_gaq.push(['_trackPageview']);
(function() {
var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
})();
//Lycos Init
function getReferrer() {
var all= this.document.cookie;
if (all== '') return false;
var cookie_name = 'REFERRER=';
var start = all.lastIndexOf(cookie_name);
if (start == -1) return false;
start += cookie_name.length;
var end = all.indexOf(';', start);
if (end == -1) end = all.length;
return all.substring(start, end);
}
function getQuery() {
var rfr = getReferrer();
if (rfr == '') return false;
var q = extractQuery(rfr, 'yahoo.com', 'p=');
if (q) return q;
q = extractQuery(rfr, '', 'q=');
return q ? q : "";
}
function extractQuery(full, site, q_param) {
var start = full.lastIndexOf(site);
if (start == -1) return false;
start = full.lastIndexOf(q_param);
if (start == -1) return false;
start += q_param.length;
var end = full.indexOf('&', start);
if (end == -1) end = full.length;
return unescape(full.substring(start, end)).split(" ").join("+");
}
function generateHref(atag, template){
atag.href=template.replace('_MYURL_', window.location.href.replace('http://', '')).replace('_MYTITLE_','Check%20out%20this%20Tripod%20Member%20site!'); 
}
var lycos_ad = Array();
var lycos_onload_timer;
var cm_role = "live";
var cm_host = "tripod.lycos.com";
var cm_taxid = "/memberembedded";
var tripod_member_name = "agentq314";
var tripod_member_page = "agentq314/Verlet/src/Verlet.cxx";
var tripod_ratings_hash = "1771502374:68d207b5df39d009f3aa6def8651af42";

var lycos_ad_category = {"dmoz":"computers\/algorithms","ontarget":"&CAT=technology&L2CAT=computing","find_what":"شرح c"};

var lycos_ad_remote_addr = "209.202.244.9";
var lycos_ad_www_server = "www.tripod.lycos.com";
var lycos_ad_track_small = "";
var lycos_ad_track_served = "";
var lycos_search_query = getQuery();
</script>

<script type="text/javascript" src="https://scripts.lycos.com/catman/init.js"></script>

<script type='text/javascript'>
 var googletag = googletag || {};
 googletag.cmd = googletag.cmd || [];
 (function() {
   var gads = document.createElement('script');
   gads.async = true;
   gads.type = 'text/javascript';
   var useSSL = 'https:' == document.location.protocol;
   gads.src = (useSSL ? 'https:' : 'http:') +
     '//www.googletagservices.com/tag/js/gpt.js';
   var node = document.getElementsByTagName('script')[0];
   node.parentNode.insertBefore(gads, node);
 })();
</script>

<script type='text/javascript'>
 googletag.cmd.push(function() {
   googletag.defineSlot('/95963596/TRI_300X250_dfp', [300, 250], 'div-gpt-ad-1450204159126-0').addService(googletag.pubads());
   googletag.defineSlot('/95963596/TRI_above_728x90_dfp', [728, 90], 'div-gpt-ad-1450204159126-1').addService(googletag.pubads());
   googletag.defineSlot('/95963596/TRI_below_728x90_dfp', [728, 90], 'div-gpt-ad-1450204159126-2').addService(googletag.pubads());
   googletag.pubads().enableSingleRequest();
   googletag.enableServices();
 });
</script>


<script type="text/javascript"> 
(function(isV)
{
    if( !isV )
    {
        return;
    }
    var adMgr = new AdManager();
    var lycos_prod_set = adMgr.chooseProductSet();
    var slots = ["leaderboard", "leaderboard2", "toolbar_image", "toolbar_text", "smallbox", "top_promo", "footer2", "slider"];
    var adCat = this.lycos_ad_category;
    adMgr.setForcedParam('page', (adCat && adCat.dmoz) ? adCat.dmoz : 'member');
    if (this.lycos_search_query)
    {
        adMgr.setForcedParam("keyword", this.lycos_search_query);
    } 
    else if(adCat && adCat.find_what)
    {
        adMgr.setForcedParam('keyword', adCat.find_what);
    }
    
    for (var s in slots)
    {
        var slot = slots[s];
        if (adMgr.isSlotAvailable(slot))
        {
            this.lycos_ad[slot] = adMgr.getSlot(slot);
        }
    }

    adMgr.renderHeader();
    adMgr.renderFooter();
}((function() {

var w = 0, h = 0, minimumThreshold = 300;

if (top == self)
{
    return true;
}
if (typeof(window.innerWidth) == 'number' )
{
    w = window.innerWidth;
    h = window.innerHeight;
}
else if (document.documentElement && (document.documentElement.clientWidth || document.documentElement.clientHeight))
{
    w = document.documentElement.clientWidth;
    h = document.documentElement.clientHeight;
}
else if (document.body && (document.body.clientWidth || document.body.clientHeight))
{
    w = document.body.clientWidth;
    h = document.body.clientHeight;
}
return ((w > minimumThreshold) && (h > minimumThreshold));
}())));




window.onload = function()
{
    var f = document.getElementById("FooterAd");
    var b = document.getElementsByTagName("body")[0];
    b.appendChild(f);
    f.style.display = "block";
    document.getElementById('lycosFooterAdiFrame').src = '/adm/ad/footerAd.iframe.html';
    


    
    // DOM Inj Ad
    (function(isTrellix)
    {
        var e = document.createElement('iframe');
        e.style.border = '0';
        e.style.margin = 0;
        e.style.display = 'block';
        e.style.cssFloat = 'right';
        e.style.height = '254px';
        e.style.overflow = 'hidden';
        e.style.padding = 0;
        e.style.width = '300px';


        var isBlokedByDomain = function( href )
        {
            var blockedDomains = [
                "ananyaporn13000.tripod.com",
                "xxxpornxxx.tripod.com"
            ];
            var flag = false;
            
            for( var i=0; i<blockedDomains.length; i++ )
            {
                if( href.search( blockedDomains[ i ] ) >= 0 )
                {
                    flag = true;
                }
            }
            return flag;
        }

        var getMetaContent = function( metaName )
        {
            var metas = document.getElementsByTagName('meta');
            for (i=0; i<metas.length; i++)
            { 
                if( metas[i].getAttribute("name") == metaName )
                { 
                    return metas[i].getAttribute("content"); 
                } 
            }
            return false;
        }
        
        var getCommentNodes = function(regexPattern)
        {
            var nodes = {};
            var nodesA = [];
            var preferredNodesList = ['a', 'c', 'b'];
        
            (function getNodesThatHaveComments(n, pattern)
            {
                if (n.hasChildNodes())
                {
                    if (n.tagName === 'IFRAME')
                    {
                        return false;
                    }
                    for (var i = 0; i < n.childNodes.length; i++)
                    {
                        if ((n.childNodes[i].nodeType === 8) && (pattern.test(n.childNodes[i].nodeValue)))
                        {
                            var areaName = pattern.exec(n.childNodes[i].nodeValue)[1];
                            nodes[areaName] = n;
                        }
                        else if (n.childNodes[i].nodeType === 1)
                        {
                            getNodesThatHaveComments(n.childNodes[i], pattern);
                        }
                    }
                }
            }(document.body, regexPattern));

            for (var i in preferredNodesList)
            {
                if (nodes[preferredNodesList[i]])
                {
                    if( isTrellix && nodes[preferredNodesList[i]].parentNode.parentNode.parentNode.parentNode )
                    {
                        nodesA.push(nodes[preferredNodesList[i]].parentNode.parentNode.parentNode.parentNode);
                    }
                    else
                    {
                        nodesA.push( nodes[preferredNodesList[i]] );
                    }
                }
            }
            return nodesA;
        }
        
        
        var properNode = null;
        var areaNodes = getCommentNodes( new RegExp( '^area Type="area_(\\w+)"' ) );

        for (var i = 0; i < areaNodes.length; i++)
        {
            var a = parseInt(getComputedStyle(areaNodes[i]).width);
            if ((a >= 300) && (a <= 400))
            {
                properNode = areaNodes[i];
                break;
            }
        }


        var propertyName = getMetaContent("property") || false;
        if( isTrellix && (properNode) )
        {
            e.src = '/adm/ad/injectAd.iframe.html';
            properNode.insertBefore(e, properNode.firstChild);
        }
        else if( isTrellix && !( properNode ) ) // Slap the ad eventhought there is no alocated slot
        {
            e.src = '/adm/ad/injectAd.iframe.html';
            e.style.cssFloat = 'none';
            var cdiv = document.createElement('div');
            cdiv.style = "width:300px;margin:10px auto;";
            cdiv.appendChild( e );
            b.insertBefore(cdiv, b.lastChild);
        }
        else if( !isBlokedByDomain( location.href ) )
        {
            var injF = document.createElement('iframe');
            injF.style.border = '0';
            injF.style.margin = 0;
            injF.style.display = 'block';
            injF.style.cssFloat = 'none';
            injF.style.height = '254px';
            injF.style.overflow = 'hidden';
            injF.style.padding = 0;
            injF.style.width = '300px';
            injF.src = '/adm/ad/injectAd.iframe.html';

            if( b && ( !isTrellix || ( typeof isTrellix == "undefined" ) ) ) // All other tripod props
            {
                var cdiv = document.createElement('div');
                cdiv.style = "width:300px;margin:10px auto;";
                cdiv.appendChild( injF );
                b.insertBefore(cdiv, b.lastChild);
            } 
        }
  }( document.isTrellix ));
}

</script>

<div id="tb_container" style="background:#DFDCCF; border-bottom:1px solid #393939; position:relative; z-index:999999999!important">
<!--form name="search" onSubmit="return searchit()" id='header_search' >
<input type="text" placeholder="Search" size=30 name="search2" value="">
<input type="button" value="Go!" onClick="searchit()">
</form>
<style>
form#header_search {
    width: 916px;
    margin: 0 auto 8px;
    position: relative;
}


form#header_search input {
    height: 40px;
    font-size: 14px;
    line-height: 40px;
    padding: 0 8px;
    box-sizing: border-box;
    background: #F4F2E9;
    border: 1px solid #BBB8B8;
    transition: background-color 300ms ease-out,
                color 300ms ease;
}

form#header_search input[type="text"] {
    width: 100%;
}
form#header_search input[type="text"]:focus {
    border-color: #A2D054;
    background-color: #fff;
    box-shadow: 0 0px 12px -4px #A2D054;
}



form#header_search input[type="button"] {
    position: absolute;
    top: 1px;
    right: 1px;
    opacity: 1;
    background: #DFDCCF;
    color: #463734;
    width: 125px;
    cursor: pointer;
    height: 38px;
    border: none;
}
form#header_search input[type="text"]:focus ~ input[type='button']:hover,
form#header_search input[type='button']:hover {
    background-color: #A5CE56;
    color: #fff;
}
form#header_search input[type="text"]:focus ~ input[type='button'] {
    background-color: #52AEDF;
    color: #fff;
}

</style>

<script>
function searchit(){
    
    // determine environment 
    var search_env 
    if (lycos_ad_www_server.indexOf(".pd.") > -1) {
       search_env = 'http://search51.pd.lycos.com/a/';
    } else if (lycos_ad_www_server.indexOf(".qa.") > -1) {
       search_env = 'http://search51.qa.lycos.com/a/';
    } else {
       search_env = 'http://search51.lycos.com/a/';
    }

var search_term = encodeURIComponent(document.search.search2.value)
var search_url = search_env+search_term;
window.open(search_url);

return false
}
</script-->
<style>
    .adCenterClass{margin:0 auto}
</style>
<div id="tb_ad" class="adCenterClass" style="display:block!important; overflow:hidden; width:916px;">

<div id="ad_container" style="display:block!important; float:left; width:728px ">
<script type="text/javascript">
if (typeof lycos_ad !== "undefined" && "leaderboard" in lycos_ad) {
  document.write(lycos_ad['leaderboard']);
}
</script>
</div>
</div>
</div>
<script type="text/javascript">
if (typeof lycos_ad !== "undefined" && "slider" in lycos_ad) {
  document.write(lycos_ad['slider']);
}
</script> <!-- added 7/22 -->
<div id="FooterAd" style="background:#DFDCCF; border-top:1px solid #393939; clear:both; display:none; width:100%!important; position:relative; z-index:999999!important; height:90px!important"> 
<div class="adCenterClass" style="display:block!important; overflow:hidden; width:916px;">
<div id="footerAd_container" style="display:block!important; float:left; width:728px">
<iframe id="lycosFooterAdiFrame" style="border:0; display:block; float:left; height:96px; overflow:hidden; padding:0; width:750px"></iframe>
</div>
</div>
</div>