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!