Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran quassiclassical trajectory help

  1. Jun 11, 2012 #1
    Hello All,

    This is my first post on any forum so I will try not to get too fancy. I have been writing in fortran for about two months now and it is my first programming experience. First, I wrote a classical dynamics program for a particle moving toward a potential barrier. Then, I wrote a dynamics program to model a harmonic oscillator with a constant frequency. Then I put the two together (but uncoupled) so that there is a harmonic oscillator moving as well as a particle approaching the potential barrier in the same simulation.
    Now I am trying to write a program that couples the oscillator to the motion of the particle. So basically it would be a diatomic molecule approaching some potential barrier and if it has enough energy, products will be formed (change in oscillator frequency). If there isn't enough energy in the initial conditions, the diatomic will bounce off the barrier back to reactants.

    The problem: I am having trouble getting energy to be conserved in the coupled program and I think it stems in my description of the oscillator and its frequency or its derivative, but I've been through it and am stumped.

    With all that said, the way I couple the oscillator to the reaction path is through a switching function, called w_q, which is coupled by transferring the oscillator motion into translational motion along the reaction path. This function is also in charge of determining the frequency of the oscillator. If reactants, the frequency is one number, and if the diatomic reacts to products, the vibrational frequency will switch to another number. All of this is dependent on the reaction path, s.

    To define some terms: s = reaction path, p_s = momentum of particle,
    q = oscillator displacement, p_q = oscillator momentum, H = total energy

    The integration method is a 4th order Runge-Kutta and the subroutine is just the derivatives to integrate newton's equations of motion.

    My input file (cpld1_in):

    1 !ntraj
    100 !nprint
    -6.01 !sreflect boundary condition
    4.0 !sreact boundary condition
    2 !nee number of energies to propagate
    1.0 !alpha
    1.0 !A barrier height
    16.0 !m molecule mass
    1.0 !mu reduced mass
    0.8 !emin minimum energy attempt
    1.10 !emax maximum energy attempt
    -6.0 !sinit reactant position
    2.0 !q_max maximum displacement
    0.005 !dt

    My fortran code (cpld20.f):
    -------------------------------------------

    program coupled20

    implicit none

    common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev

    real(8) :: p_q, pq_2, pq_3, pq_4
    real(8) :: p_q_dot, dpqdt2, dpqdt3, dpqdt4
    real(8) :: q, q_2, q_3, q_4, q_dot, dqdt2, dqdt3, dqdt4
    real(8) :: alpha, A, s_0, m, mu, dt, tutoev
    real(8) :: rand, sinit, q_max, sref, sreact
    real(8) :: einit, emin, emax, V, K, H, w_q, phi, pi, hbar
    real(8) :: sk1, sk2, sk3, sk4, psk1, psk2, psk3, psk4
    real(8) :: qk1, qk2, qk3, qk4, pqk1, pqk2, pqk3, pqk4
    real(8) :: s, s_2, s_3, s_4, s_dot, dsdt2, dsdt3, dsdt4
    real(8) :: p_s,ps_2,ps_3,ps_4,p_s_dot,dpsdt2,dpsdt3,dpsdt4

    integer :: i,ntraj,nen,nee,nstep,nprint

    open(1,file='cpld1_in',status='unknown')

    read(1,*)ntraj
    read(1,*)nprint
    read(1,*)sref
    read(1,*)sreact
    read(1,*)nee
    read(1,*)alpha
    read(1,*)A
    read(1,*)m
    read(1,*)mu
    read(1,*)emin
    read(1,*)emax
    read(1,*)sinit
    read(1,*)q_max
    read(1,*)dt

    tutoev = 1.0364314
    pi = acos(-1e0)
    s_0 = 0.0


    open(unit=2, file='cpld20_out', action='write', status='replace')
    write(2,*) 'Time(fs) s(A) Displace(A) PE(eV) KE(eV) H(eV)'

    do nen=1,nee

    einit=emin+(emax-emin)*float(nen-1)/float(nee-1)
    call random_number(rand)
    phi = 2.0*pi*rand

    do i=1,ntraj

    s = sinit

    w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18

    p_s = sqrt(2.0*m*einit/tutoev)
    q = q_max*cos(phi)
    p_q = -mu*w_q*q_max*sin(phi)

    nstep = 0

    do while(s > sref .and. s < sreact)

    if(mod(nstep,nprint) == 0) then
    V = A*exp(-alpha*(s-s_0)**2) + mu*w_q*w_q*q*q/2.0
    K = tutoev*p_s*p_s/(2.0*m) + p_q*p_q/(2.0*mu)
    H = V + K
    write(2,130)10.0*dt*nstep,s,q,V,K,H
    end if

    call cpld_der(s,p_s,q,p_q,s_dot,p_s_dot,q_dot,p_q_dot)

    sk1 = dt*s_dot
    psk1 = dt*p_s_dot
    qk1 = dt*q_dot
    pqk1 = dt*p_q_dot
    s_2 = s + sk1/2.0
    ps_2 = p_s + psk1/2.0
    q_2 = q + qk1/2.0
    pq_2 = p_q + pqk1/2.0

    call cpld_der(s_2,ps_2,q_2,pq_2,dsdt2,dpsdt2,dqdt2,dpqdt2)

    sk2 = dt*dsdt2
    psk2 = dt*dpsdt2
    qk2 = dt*dqdt2
    pqk2 = dt*dpqdt2
    s_3 = s + sk2/2.0
    ps_3 = p_s + psk2/2.0
    q_3 = q + qk2/2.0
    pq_3 = p_q + pqk2/2.0

    call cpld_der(s_3,ps_3,q_3,pq_3,dsdt3,dpsdt3,dqdt3,dpqdt3)

    sk3 = dt*dsdt3
    psk3 = dt*dpsdt3
    qk3 = dt*dqdt3
    pqk3 = dt*dpqdt3
    s_4 = s + sk3
    ps_4 = p_s + psk3
    q_4 = q + qk3
    pq_4 = p_q + pqk3

    call cpld_der(s_4,ps_4,q_4,pq_4,dsdt4,dpsdt4,dqdt4,dpqdt4)

    sk4 = dt*dsdt4
    psk4 = dt*dpsdt4
    qk4 = dt*dqdt4
    pqk4 = dt*dpqdt4
    s = s + (sk1 + 2.0*sk2 + 2.0*sk3 + sk4)/6.0
    p_s = p_s + (psk1 + 2.0*psk2 + 2.0*psk3 + psk4)/6.0
    q = q + (qk1 + 2.0*qk2 + 2.0*qk3 + qk4)/6.0
    p_q = p_q + (pqk1 + 2.0*pqk2 + 2.0*pqk3 + pqk4)/6.0

    w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18

    nstep = nstep + 1
    end do
    end do
    end do

    130 format(6f9.5)
    close(2)

    end program coupled20

    subroutine cpld_der(s,ps,q,pq,sdot,psdot,qdot,pqdot)

    implicit none

    common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev

    real(8) :: A, alpha, m, mu, w_q, q_max, phi, s_0, tutoev
    real(8) :: s, ps, q, pq, sdot, psdot, qdot, pqdot
    real(8) :: dwds1,dwds2,dw2ds,psdt1,psdt2,psdt3

    sdot = ps/m

    dwds1 = -0.19*alpha*exp(alpha*(s-s_0))
    dwds2 = dwds1/(1.0+exp(alpha*(s-s_0)))**2
    dw2ds = 2.0*dwds2*w_q

    psdt1 = -pq*q_max*sin(phi)*dwds2
    psdt2 = -2.0*alpha*(s-s_0)*A*exp(-alpha*(s-s_0)**2)
    psdt3 = 0.5*mu*q*q*dw2ds
    psdot = -1.0*(psdt1+psdt2+psdt3)

    qdot = pq/mu
    pqdot = -mu*w_q**2*q

    end subroutine cpld_der
    ------------------------------------------------------------------

    Questions are welcome and I really appreciate the effort of anyone who has read this far.
    If I should reformat my post, let me know.
    Thank you!
     
  2. jcsd
  3. Jun 12, 2012 #2
    I'm not going to go through line by line, but here are some recommendations:

    1. make more meaningful variable names (reac_pth versus s makes this easier to read; do not be afraid of continuations)
    2. make a data dictionary where you explain your variables so the equations are more clear
    3. insert sanity checks. Toss in a write statement coupled with a pause so that you see some output, check it against what you know should be happening, and continue. This will allow you to see where the issue is.

    Cheers
     
  4. Jun 12, 2012 #3
    Thanks for your help swartzism.
     
  5. Jun 12, 2012 #4
    So now that I have a little more time,

    Code (Text):
    program coupled20

    implicit none

    common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev

    real(8) :: p_q, pq_2, pq_3, pq_4
    real(8) :: p_q_dot, dpqdt2, dpqdt3, dpqdt4
    real(8) :: q, q_2, q_3, q_4, q_dot, dqdt2, dqdt3, dqdt4
    real(8) :: alpha, A, s_0, m, mu, dt, tutoev
    real(8) :: rand, sinit, q_max, sref, sreact
    real(8) :: einit, emin, emax, V, K, H, w_q, phi, pi, hbar
    real(8) :: sk1, sk2, sk3, sk4, psk1, psk2, psk3, psk4
    real(8) :: qk1, qk2, qk3, qk4, pqk1, pqk2, pqk3, pqk4
    real(8) :: s, s_2, s_3, s_4, s_dot, dsdt2, dsdt3, dsdt4
    real(8) :: p_s,ps_2,ps_3,ps_4,p_s_dot,dpsdt2,dpsdt3,dpsdt4

    integer :: i,ntraj,nen,nee,nstep,nprint

    open(1,file='cpld1_in',status='unknown')

    read(1,*)ntraj
    read(1,*)nprint
    read(1,*)sref
    read(1,*)sreact
    read(1,*)nee
    read(1,*)alpha
    read(1,*)A
    read(1,*)m
    read(1,*)mu
    read(1,*)emin
    read(1,*)emax
    read(1,*)sinit
    read(1,*)q_max
    read(1,*)dt

    [B]print*, <above vars> ! This will allow you to see if everything is read in as you want[/B]
    [B]pause[/B]

    tutoev = 1.0364314
    pi = acos(-1e0)
    s_0 = 0.0


    open(unit=2, file='cpld20_out', action='write', status='replace')
    write(2,*) 'Time(fs) s(A) Displace(A) PE(eV) KE(eV) H(eV)'

    [B]! Here, I suggest labeling these loops with do.. continue vs do.. enddo; also, indent[/B]
    do 100 nen=1,nee

       einit=emin+(emax-emin)*float(nen-1)/float(nee-1)
       call random_number(rand)
       phi = 2.0*pi*rand

       [B]print*, <above vals>[/B]
       [B]pause[/B]

       do 200 i=1,ntraj

          s = sinit

          w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18

          p_s = sqrt(2.0*m*einit/tutoev)
          q = q_max*cos(phi)
          p_q = -mu*w_q*q_max*sin(phi)
         
          [B]print*, <above vals>[/B]
          [B]pause[/B]

          nstep = 0

          do 300 while(s > sref .and. s < sreact)

             if(mod(nstep,nprint) == 0) then
                V = A*exp(-alpha*(s-s_0)**2) + mu*w_q*w_q*q*q/2.0
                K = tutoev*p_s*p_s/(2.0*m) + p_q*p_q/(2.0*mu)
                H = V + K
                write(2,130)10.0*dt*nstep,s,q,V,K,H
             end if

             call cpld_der(s,p_s,q,p_q,s_dot,p_s_dot,q_dot,p_q_dot)

             sk1 = dt*s_dot
             psk1 = dt*p_s_dot
             qk1 = dt*q_dot
             pqk1 = dt*p_q_dot
             s_2 = s + sk1/2.0
             ps_2 = p_s + psk1/2.0
             q_2 = q + qk1/2.0
             pq_2 = p_q + pqk1/2.0

             [B]print*, <above vals>[/B]
             [B]pause[/B]

             call cpld_der(s_2,ps_2,q_2,pq_2,dsdt2,dpsdt2,dqdt2,dpqd t2)

             sk2 = dt*dsdt2
             psk2 = dt*dpsdt2
             qk2 = dt*dqdt2
             pqk2 = dt*dpqdt2
             s_3 = s + sk2/2.0
             ps_3 = p_s + psk2/2.0
             q_3 = q + qk2/2.0
             pq_3 = p_q + pqk2/2.0

             [B]print*, <above vals>[/B]
             [B]pause[/B]

             call cpld_der(s_3,ps_3,q_3,pq_3,dsdt3,dpsdt3,dqdt3,dpqd t3)

             sk3 = dt*dsdt3
             psk3 = dt*dpsdt3
             qk3 = dt*dqdt3
             pqk3 = dt*dpqdt3
             s_4 = s + sk3
             ps_4 = p_s + psk3
             q_4 = q + qk3
             pq_4 = p_q + pqk3

             [B]print*, <above vals>[/B]
             [B]pause[/B]

             call cpld_der(s_4,ps_4,q_4,pq_4,dsdt4,dpsdt4,dqdt4,dpqd t4)

    [B]! Tired of spacing in.. my preference is 3 per new indentation[/B]

    sk4 = dt*dsdt4
    psk4 = dt*dpsdt4
    qk4 = dt*dqdt4
    pqk4 = dt*dpqdt4
    s = s + (sk1 + 2.0*sk2 + 2.0*sk3 + sk4)/6.0
    p_s = p_s + (psk1 + 2.0*psk2 + 2.0*psk3 + psk4)/6.0
    q = q + (qk1 + 2.0*qk2 + 2.0*qk3 + qk4)/6.0
    p_q = p_q + (pqk1 + 2.0*pqk2 + 2.0*pqk3 + pqk4)/6.0

    w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18

             [B]print*, <above vals>[/B]
             [B]pause[/B]

    nstep = nstep + 1
    [B]
     300       continue
     200    continue
     100 continue[/B]

    130 format(6f9.5)
    close(2)

    end program coupled20

    subroutine cpld_der(s,ps,q,pq,sdot,psdot,qdot,pqdot)

    implicit none

    common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev

    real(8) :: A, alpha, m, mu, w_q, q_max, phi, s_0, tutoev
    real(8) :: s, ps, q, pq, sdot, psdot, qdot, pqdot
    real(8) :: dwds1,dwds2,dw2ds,psdt1,psdt2,psdt3

    sdot = ps/m

    dwds1 = -0.19*alpha*exp(alpha*(s-s_0))
    dwds2 = dwds1/(1.0+exp(alpha*(s-s_0)))**2
    dw2ds = 2.0*dwds2*w_q

    psdt1 = -pq*q_max*sin(phi)*dwds2
    psdt2 = -2.0*alpha*(s-s_0)*A*exp(-alpha*(s-s_0)**2)
    psdt3 = 0.5*mu*q*q*dw2ds
    psdot = -1.0*(psdt1+psdt2+psdt3)

    qdot = pq/mu
    pqdot = -mu*w_q**2*q

    end subroutine cpld_der
    So keep in mind, none of this is exact. Put the print & pause combos in places you actually know the values (sanity checks), format a little so its easier to read (my suggestion is the do.. continue). Also, try inserting the equations in full form via comment above each instance so the reader knows where you are in your calculations and can easily see if there is something fishy in your math. Including references helps.

    I compile with gfortran; your compiler might tell you that it deleted the pauses. Fear not, it will still pause and ask for you to press enter or type "go" to continue. Essentially a debugging break point.
     
  6. Jun 12, 2012 #5
    My apologies for the lack of comments and indents. When I pasted the code into my post above, it took out my indents.? I do like a clean code in that sense.
    When it comes to the comments, I had a few essentials in there, but after about the 10th variation of the code, I started leaving out the comments and then when I wrote my post it slipped my mind to add them back in.

    Thanks for showing me the do.. continue. It is a nice label for easy reading as opposed to the end do style. I didn't know you could add pauses in like that, which helps me see how printing output to the display make a lot more sense in some cases than writing to a file.

    Most of today was spent doing checks on the physics by hand, so it would be a big help incorporating the sanity checks into the code.
    I checked all the physics on my previous versions with the uncoupled oscillator and particle and found that the initial conditions agree with my coupled version, which is good considering that all the initial parameters are the same.
    This leads me to believe that the error is in the derivative somewhere because I used this exact same integration method on the uncoupled version, which conserved energy.
    Either way, all of the cleaning up and sanity checks will help.

    After I clean it up I may repost.

    Thanks again!
     
  7. Jun 12, 2012 #6

    jtbell

    User Avatar

    Staff: Mentor

    To preserve indentation, enclose your code in [noparse]
    Code (Text):
     
    [/noparse] tags.
     
  8. Jun 26, 2012 #7
    Thanks everyone for your help! Found an error in one of my derivatives for the equations of motion and all is well. You advice was still used as I created a dictionary for terms, added more comments and optimized the outputs to help with finding errors.
    I'm not sure how this works, but any moderators reading this may close this thread at their own will.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Fortran quassiclassical trajectory help
  1. Fortran help (Replies: 0)

  2. Help in fortran (Replies: 2)

  3. [Fortran] help (Replies: 2)

  4. FORTRAN Help (Replies: 1)

  5. Fortran Help (Replies: 2)

Loading...