Fortran quassiclassical trajectory help

[movezig]

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!

 

[swartzism]

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

 

[movezig]

Thanks for your help swartzism.

 

[swartzism]

So now that I have a little more time,

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.0open(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.

 

[movezig]

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!

 

[jtbell]

When I pasted the code into my post above, it took out my indents.?

To preserve indentation, enclose your code in [noparse]

[/noparse] tags.

 

[movezig]

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.