One Dimensional System of the LJ Model

If you would like to work on your LJ Model only on one dimension, you will need to modify this part of the program. Also your input file should be modified. See demonstration below:

Click here to see example of Input file

Before:

Molec_dyn.f
* Read initial positions, velocities
      DO  I = 1,NATOM
         READ(2,*) X(I), Y(I), Z(I)
         READ(2,*) U(I), V(I), W(I)
      END DO

  DO I = 1,NATOM
         KE = KE + (U(I)**2 + V(I)**2 + W(I)**2)/2
      END DO
      WRITE(3,99) K*0.0, 0, PE, KE, PE+KE
      DO I = 1,NATOM
         WRITE(3,99) K*H, I, X(I), Y(I), Z(I)
         WRITE(3,99) K*H, I, U(I), V(I), W(I)
      END DO

* Main loop
      DO K = 1, NSTEP
         DO I = 1,NATOM
            X(I) = X(I) + H*(U(I) + HO2*FX(I,B1))
            Y(I) = Y(I) + H*(V(I) + HO2*FY(I,B1))
            Z(I) = Z(I) + H*(W(I) + HO2*FZ(I,B1))
         END DO
         CALL FORCE_LJ(NATOM,NATOMX,X,Y,Z,B2,FX,FY,FZ)
         DO I = 1,NATOM
            U(I) = U(I) + HO2*(FX(I,B1) + FX(I,B2))
            V(I) = V(I) + HO2*(FY(I,B1) + FY(I,B2))
            W(I) = W(I) + HO2*(FZ(I,B1) + FZ(I,B2))
         END DO

Potentional_lj.f
         IF (MOD(K,PRNFREQ) .EQ. 0) THEN
            PE = POTENTIAL_LJ(NATOM,X,Y,Z)
            KE = 0
            DO I = 1,NATOM
               KE = KE + (U(I)**2 + V(I)**2 + W(I)**2)/2
            END DO
            WRITE(3,99)INT(1000.*K*H),0,PE,KE,PE+KE
            DO I = 1,NATOM
               WRITE(3,99)K*H,I,X(I),Y(I),Z(I)
               WRITE(3,99)K*H,I,U(I),V(I),W(I)
            END DO
         END IF

* Compute potential.
      DO I = 1, NATOM-1
         DO J = I+1, NATOM
            RIJSQ = 1.0/((X(I)-X(J))**2 + (Y(I)-Y(J))**2 +
     $           (Z(I)-Z(J))**2)
            RIJM6 = RIJSQ*RIJSQ*RIJSQ
            POT = POT + RIJM6*(RIJM6 - 2.0)
                    END DO
      END DO

After:

Molec_dyn.f
* Read initial positions, velocities
      DO  I = 1,NATOM
         READ(2,*) X(I)
         READ(2,*) U(I)
      END DO

  DO I = 1,NATOM
         KE = KE + (U(I)**2 )/2
      END DO
      WRITE(3,99) K*0.0, 0, PE, KE, PE+KE
      DO I = 1,NATOM
         WRITE(3,99) K*H, I, X(I)
         WRITE(3,99) K*H, I, U(I)
      END DO

* Main loop
      DO K = 1, NSTEP
         DO I = 1,NATOM
            X(I) = X(I) + H*(U(I) + HO2*FX(I,B1))
         END DO
         CALL FORCE_LJ(NATOM,NATOMX,X,Y,Z,B2,FX,FY,FZ)
         DO I = 1,NATOM
            U(I) = U(I) + HO2*(FX(I,B1) + FX(I,B2))
         END DO

Potentional_lj.f
         IF (MOD(K,PRNFREQ) .EQ. 0) THEN
            PE = POTENTIAL_LJ(NATOM,X,Y,Z)
            KE = 0
            DO I = 1,NATOM
               KE = KE + (U(I)**2)/2
            END DO
            WRITE(3,99)INT(1000.*K*H),0,PE,KE,PE+KE
            DO I = 1,NATOM
               WRITE(3,99)K*H,I,X(I)
               WRITE(3,99)K*H,I,U(I)
            END DO
         END IF

* Compute potential.
      DO I = 1, NATOM-1
         DO J = I+1, NATOM
            RIJSQ = 1.0/((X(I)-X(J))**2 
                RIJM6 = RIJSQ*RIJSQ*RIJSQ
            POT = POT + RIJM6*(RIJM6 - 2.0)
                    END DO
      END DO