anharmonic.xmds

<?xml version="1.0"?>
<!--Example damped anharmonic oscilator simulation using nu-P (+P with squeezed-state basis) function-->
<!--Essentially a complicated set of coupled stochastic ODEs-->
<!--Noise correlations are especially complicated, which slows things down-->
<!--The equations are also not very stable (no gauges)-->

<simulation>

  <name>anharmonic</name>

  <author>Joel Corney</author>
  <description>
    Example damped anharmonic oscilator simulation using nu-P (+P with
    squeezed-state basis) function.  Essentially a complicated set of
    coupled stochastic ODEs.  Noise correlations are especially
    complicated, which slows things down.  The equations are also not
    very stable (no gauges)
  </description>

  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>1000</paths>
  <use_mpi>yes</use_mpi>
  <seed>1 2</seed>
  <noises>14</noises>


  /***********/
  <globals>
  /***********/
  <![CDATA[
    const double Pi = 2.0*asin(1.0);
    const double N = 1000;         /*Number of particles*/
    const double kappa = 0.5/N;    /*Interaction strength*/
    const double gamm = 100*kappa; /*Damping rate*/
  ]]>
  </globals>

  /*************/
  <field>
  /*************/
    <name>main</name>
    <samples>1</samples>
    
    /***MAIN***/
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>alpha1 alpha2 nu1 nu2 mu</components>
      <fourier_space>no</fourier_space>
      
      <![CDATA[
        /* NB: We take the '2' variables to be the conjugates */
        alpha1 = rcomplex(sqrt(N),0.0);
        alpha2 = rcomplex(sqrt(N),0.0);    
        nu1 = rcomplex(0.0,0.0);
        nu2 = rcomplex(0.0,0.0);
        mu = c_sqrt(1+ nu1+nu2);
      ]]>
    </vector>        
  </field>

  /******************/
  <sequence>
  /******************/
 
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>5</interval>
      <lattice>5000</lattice>
      <samples>200</samples>
      <iterations>4</iterations>
      <vectors>main</vectors>
      <![CDATA[
        complex nunu = nu1*nu2;

        complex w5 = rcomplex(n_5,n_6)/sqrt(2.0);       /*Complex noise*/
        complex w6 = rcomplex(n_7,n_8)/sqrt(2.0);       /*Complex noise*/
        complex w7 = rcomplex(n_9,n_10)/sqrt(2.0);      /*Complex noise*/
        complex w8 = rcomplex(n_11,n_12)/sqrt(2.0);     /*Complex noise*/
        complex w9 = rcomplex(n_13,n_14)/sqrt(2.0);     /*Complex noise*/

        complex s1 = c_sqrt(8*i*kappa*alpha1*alpha2*nu1*mu + gamm*mu*nu1);      /*Noise Correlations*/
        complex s2 = c_sqrt(-8*i*kappa*alpha1*alpha2*nu2*mu + gamm*mu*nu2);     /*Noise Correlations*/
        complex s3 = c_sqrt(-8*i*kappa*nu1*nu1*mu*mu);  /*Noise Correlations*/
        complex s4 = c_sqrt(8*i*kappa*nu2*nu2*mu*mu);   /*Noise Correlations*/

        complex s5 = c_sqrt(2*i*kappa*nu1*(2*alpha2*nu1*mu - alpha1*(3*nu1*nu2 + 2)));  /*Noise Correlations*/
        complex s6 = c_sqrt(2*i*kappa*nu1*nu2*(alpha1*nu2 - 2*alpha2*mu));  /*Noise Correlations*/
        complex s7 = c_sqrt(-2*i*kappa*nu1*nu2*(alpha2*nu1 - 2*alpha1*mu));  /*Noise Correlations*/
        complex s8 = c_sqrt(-2*i*kappa*nu2*(2*alpha1*nu2*mu - alpha2*(3*nu1*nu2 + 2)));  /*Noise Correlations*/
        complex s9 = c_sqrt(-gamm*nu1*nu2);

        complex nu1_deriv = i*kappa*((alpha2*alpha2*nu1 + alpha1*alpha1*(2+3*nu1*nu2))/mu
                + 3*nu1) + s3*n_3 + s5*~w5 + s7*~w7 + s9*w9;
        complex nu2_deriv = -i*kappa*((alpha1*alpha1*nu2 + alpha2*alpha2*(2+3*nu1*nu2))/mu
                + 3*nu2) + s4*n_4 + s6*~w6 + s8*~w8;

        dalpha1_dt =-2*i*kappa*(alpha1*(alpha1*alpha2  + nu1*nu2) 
                - alpha1/2 + alpha2*nu1*mu/2)- gamm*alpha1/2
                   + s1*n_1 + s5*w5 + s6*w6+ s9*w9;
        dalpha2_dt = 2*i*kappa*(alpha2*(alpha1*alpha2 + nu1*nu2) 
                - alpha2/2 + alpha1*nu2*mu/2)- gamm*alpha2/2
                   + s2*n_2 + s7*w7 + s8*w8+ s9*~w9;
        dnu1_dt = nu1_deriv;
        dnu2_dt = nu2_deriv;
        dmu_dt = (nu2*nu1_deriv + nu1*nu2_deriv)/(2*mu);
      ]]>
    </integrate>
  </sequence>

  /**************/
  <output>
    <filename>anharmonic.xsil</filename>
    <group>
      <sampling>
       <vectors>main</vectors>
        <moments>mumu nunu nunui N_alpha N_alphai alpha_x alpha_y alpha_x2</moments>
        <![CDATA[
          mumu = mu*mu;
          nunu = nu1*nu2;
          N_alpha = alpha1*alpha2 + nunu;
          nunui = imag(nu1*nu2);
          N_alphai = imag(alpha1*alpha2 +nu1*nu2);
          alpha_x = (alpha1 + alpha2)/2;
          alpha_y = (alpha1 - alpha2)/(2*i);
          alpha_x2 = (alpha1 + alpha2)*(alpha1 + alpha2)/4.0 + (1.0 + 2*nunu - mu*(nu1+nu2))/4.0;
        ]]>
      </sampling>
    </group>
   </output>
</simulation>