0

I am having some trouble with two of my modules for a molecular dynamics program I am making. both occur at the force calculation in the module and I have tried changing the precision of all of my values. One of the modules is printed below. I get the error at line 14, the first force calculation and I assume that the rest of them will also have the issue. Any suggestions?

module morseforce

contains

  subroutine morse(x1,x2,y1,y2,z1,z2,axi,axj,ayi,ayj,azi,azj,m)
    double precision,intent(in)::x1,x2,y1,y2,z1,z2
    double precision,intent(out)::axi,axj,ayi,ayj,azi,azj
    double precision::rx,ry,rz,rmag
    double precision::fx2,fy2,fz2,fx1,fy1,fz1,m
    double precision::De=118.1d0,alpha=1.253d0, re=3.83d0
    rx=x2-x1
    ry=y2-y1
    rz=z2-z1
    rmag=sqrt(rx**2.0d0+ry**2.0d0+rz**2.0d0)
    fx2=-2.0*De*(1-exp(-alpha*(rx-re))*(alpha)*exp(-alpha*(rx-re)))*(rx/rmag)
    fy2=-2.0*De*(1-exp(-alpha*(ry-re))*(alpha)*exp(-alpha*(ry-re)))*(ry/rmag)
    fz2=-2.0*De*(1-exp(-alpha*(rz-re))*(alpha)*exp(-alpha*(rz-re)))*(rz/rmag)
    fx1=-fx2
    fy1=-fy2
    fz1=-fz2
    axi=fx1/m
    axj=fx2/m
    ayi=fy1/m
    ayj=fy2/m
    azi=fz1/m
    azj=fz2/m
  end subroutine

end module morseforce
Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
Megan McCracken
  • 43
  • 1
  • 1
  • 6
  • 1
    Welcome to Stack Overflow. Please, take the [tour](https://stackoverflow.com/tour) and learn [How to Ask](https://stackoverflow.com/questions/how-to-ask). Then, present us a clear description of what are you trying to do, what you got wrong and what you expected to get, along with a [Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve). Especially what kind of error (exact error message), what are the input values for the call? – albert Oct 06 '18 at 13:47
  • I will add that you should use indentation in your code (see my edit) and you should use `implicit none`. Please use that when you prepare a better code example. – Vladimir F Героям слава Oct 06 '18 at 14:11
  • @Megan, which one is line 14 (i.e., are you counting blank lines)? Don't know if it helps, but you want to integer exponents (if possible) in the line `rmag = sqrt()`. Some compiler may convert `rx**2.0d0` into `rx*rx`. Others may do `exp(2.0d0*log(rx))`. – evets Oct 06 '18 at 15:03

0 Answers0