total_energy_bigger.xmds

<?xml version="1.0"?>
<!-- Simulate evaporative cooling of an intially -->
<!-- non-condensed gas through T_c -->

<!-- This incarnation calculates the moments needed to -->
<!-- sample the total energy of the system -->

<simulation>

  <name>total_energy_bigger</name>

  <author>Tim Vaughan</author>
  <description>
    Simulate evaporative cooling of an intially
    non-condensed gas through T_c

    This incarnation calculates the moments needed to
    sample the total energy of the system
  </description>

  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <type>stochastic</type>
  <use_mpi>yes</use_mpi>
  <paths>120</paths>
  <seed>98777214</seed>
  <noises>2</noises>

  <globals>
  <![CDATA[
    #define SIN2PIX _field1[_i1]
    #define SIN2PIY _field2[_i2]
    #define SIN2PIZ _field3[_i3]
    #define SIN50PIX _field4[_i1]
    #define SIN50PIY _field5[_i2]
    #define SIN50PIZ _field6[_i3]

    double _field1[_MAIN_LATTICE_X1];
    double _field2[_MAIN_LATTICE_X2];
    double _field3[_MAIN_LATTICE_X3];
    double _field4[_MAIN_LATTICE_X1];
    double _field5[_MAIN_LATTICE_X2];
    double _field6[_MAIN_LATTICE_X3];

    const double d0 = 10e-6;
    const double t0 = 0.1422;

    const double T = 0.55541;
    const double t_max = 1.0; // time over which trap is reduced
    const double u = 7.5398e-4;

    const double Vmax = 3.5382e3;
    const double Gmax = 7.111893e2;

    const double alpha = 1.8005;

    // Constants related to initial thermal distribution:

    const double kb = 1.3806503e-23;
    const double hbar = 1.054571596e-34;
    const double m = 1.5e-25;

    const double N0 = 1.2626; // mu=-1.4e-7*kb, JC's thesis
    const double Temp = 2.4e-7;
    const double mu = -kb*Temp*log(1+1/N0);
    const double beta = 1/(kb*Temp);

    const double deltax = 1/32;
    const double deltak = 2*M_PI;
    const double cf = pow(deltak,1.5);

    // Constants defined to decrease on-the-fly CPU load:

    const complex sqrtiu = c_sqrt(i*u);
    const double deltak2 = deltak*deltak;
    const double deltak3 = deltak*deltak*deltak;
    const double deltax2 = deltax*deltax;
    const double deltax3 = deltax*deltax*deltax;
  ]]>
  </globals>

  <field>
    <dimensions>x y z</dimensions>
    <lattice>32 32 32</lattice>
    <domains>(-0.5,0.5) (-0.5,0.5) (-0.5,0.5)</domains>
    <vector>
      <name> main </name>
      <type> complex </type>
      <components>phi1 phi2</components>
      <fourier_space>no no no</fourier_space>
  </field>

  <sequence>
    <initial>
      <samples>1 1 1 1</samples>
      <fourier_space>yes yes yes</fourier_space>
      <![CDATA[
        complex rn1 = rcomplex(n1*cf,0);
        complex rn2 = rcomplex(n2*cf,0);

        double true_k2 = (kx*kx + ky*ky + kz*kz)/(d0*d0);
        double N_th = 1/(exp(beta*(hbar*hbar*true_k2/2/m - mu)) - 1);

        // Define free-field anihiliation operator amplitude:

        complex ffanihil = c_sqrt(rcomplex(N_th/2,0))*(rn1 + i*rn2);

        // This is related to the Fourier transform of the rescaled stochastic
        // fields in the following way:

        phi1 = ffanihil*rcomplex(pow(2*M_PI,-1.5),0);
        phi2 = phi1;

        // Initialize constant fields:

        SIN2PIX = pow(sin(M_PI*(_X1_MIN+_i1*_DX1)),2.0);
        SIN2PIY = pow(sin(M_PI*(_X2_MIN+_i2*_DX2)),2.0);
        SIN2PIZ = pow(sin(M_PI*(_X3_MIN+_i3*_DX3)),2.0);
        SIN50PIX = pow(sin(M_PI*(_X1_MIN+_i1*_DX1)),50.0);
        SIN50PIY = pow(sin(M_PI*(_X2_MIN+_i2*_DX2)),50.0);
        SIN50PIZ = pow(sin(M_PI*(_X3_MIN+_i3*_DX3)),50.0);
      ]]>
    </initial>

    <integrate>
      <interval>0.5554077</interval>
      <lattice>50000</lattice>
      <samples>100 100 100 100</samples>

      <k_propagation>
        <![CDATA[
          dphi1dt = -i*(kx*kx + ky*ky + kz*kz)*0.5;
          dphi2dt = -i*(kx*kx + ky*ky + kz*kz)*0.5;
        ]]>
      </k_propagation>

      <x_propagation>
        <iterations>3</iterations>
        <functions>
          <![CDATA[
            double t_eff = t;
            if (t>t_max) t_eff = t_max;

            double V = (1-alpha*t_eff)*Vmax*(SIN2PIX + SIN2PIY + SIN2PIZ);

            double G = Gmax*(SIN50PIX + SIN50PIY + SIN50PIZ);
          ]]>
        </functions>
        <![CDATA[
          dphi1dt = -i*(u*phi1*conj(phi2) + V - i*G + sqrtiu*n1)*phi1;
          dphi2dt = -i*(u*phi2*conj(phi1) + V - i*G + sqrtiu*n2)*phi2;
        ]]>
      </x_propagation>
    </integrate>
  </sequence>

  <output>
    <filename>total_energy_bigger.xsil</filename>
    <overwrite>yes</overwrite>
    <group>
      <sampling>
        <fourier_space>yes yes yes</fourier_space>
        <lattice>32 32 0</lattice>
        <moments>n_k</moments>
        <![CDATA[
          // column density:
          n_k = phi1*conj(phi2)*deltak2;
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>yes yes yes</fourier_space>
        <lattice>0 0 0</lattice>
        <moments>N_pre F_num_pre Kx_pre Ky_pre Kz_pre KE_term_pre</moments>
        <![CDATA[
          N_pre = phi1*conj(phi2);
          F_num_pre = deltak3*phi1*phi1*conj(phi2)*conj(phi2);
          Kx_pre = N_pre*kx/d0;
          Ky_pre = N_pre*ky/d0;
          Kz_pre = N_pre*kz/d0;
          KE_term_pre = hbar*hbar/2/m/d0/d0*(kx*kx + ky*ky + kz*kz)*N_pre;
        ]]>
      </sampling>
      <post_propagation>
        <fourier_space>no</fourier_space>
        <moments>N N2 F_num Kx Ky Kz KE_term COM_P_sq</moments>
        <![CDATA[
          N = N_pre;
          N2 = N*N;
          F_num = F_num_pre;
          Kx = Kx_pre;
          Ky = Ky_pre;
          Kz = Kz_pre;
          KE_term = KE_term_pre;
          COM_P_sq = hbar*hbar*(Kx*Kx+Ky*Ky+Kz*Kz) + 2*m*KE_term;
        ]]>
      </post_propagation>
    </group>
    <group>
      <sampling>
        <fourier_space>no no no</fourier_space>
        <lattice>32 32 0</lattice>
        <moments>n_x</moments>
        <![CDATA[
          // column density in position space:
          n_x = phi1*conj(phi2)*deltak2;
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>no no no</fourier_space>
        <lattice>0 0 0</lattice>
        <moments>N_x PE_term_pre Int_term_pre X_pre Y_pre Z_pre mod_X_sq_pre</moments>
        <![CDATA[
          complex N_x = phi1*conj(phi2);

          double t_eff = t;
          if (t>t_max) t_eff = t_max;

          complex V = (1-alpha*t_eff)*hbar/t0*Vmax
          *(SIN2PIX + SIN2PIY + SIN2PIZ)*(1.0+0*i);

          PE_term_pre = V*N_x;
          Int_term_pre = 0.5*(hbar/t0*u)*N_x*N_x;
          X_pre = N_x*x*d0;
          Y_pre = N_x*y*d0;
          Z_pre = N_x*z*d0;
          mod_X_sq_pre = (x*x + y*y + z*z)*N_x*d0*d0;
        ]]>
      </sampling>
      <post_propagation>
        <fourier_space>no</fourier_space>
        <moments>PE_term Int_term X Y Z mod_X_sq NX_sq X_sq</moments>
        <![CDATA[
          complex N_sq = N_x*N_x;
          PE_term = PE_term_pre;
          Int_term = Int_term_pre;
          X = X_pre; 
          Y = Y_pre; 
          Z = Z_pre; 
          mod_X_sq = mod_X_sq_pre;
          NX_sq = X*X + Y+Y + Z*Z + mod_X_sq;
          X_sq = NX_sq/N_sq;
        ]]>
      </post_propagation>
    </group>
  </output>
</simulation>