PROGRAM SQUARE
      IMPLICIT REAL*8(A-H, O-Z)
      PARAMETER(NDIM  =50)
      PARAMETER(NDIM2 = NDIM*NDIM)
      PARAMETER(PI    =3.1415926D0)
      PARAMETER(SML   = 0.1D-5)
      DIMENSION FF(0:NDIM, 0:NDIM)
      DIMENSION XX(0:NDIM),  YY(0:NDIM)
      DIMENSION RX(0:NDIM2), RY(0:NDIM2)
      DIMENSION PX(0:NDIM, 0:NDIM), PY(0:NDIM, 0:NDIM)
C
C     ***** INPUT ATOM POSITION WITHIN A CELL *****
C
      READ(5,*)     NLAT, AA, NG
      WRITE(6,2010) NLAT
 2010 FORMAT('# LATTICE SIZE=',I4)
      WRITE(6,2011) AA
 2011 FORMAT('# LATTICE CONSTANT=',F7.4) 
      WRITE(6,2012) NG
 2012 FORMAT('# CALCULATED K-RANGE =', I4, ' x G') 
      READ(5,*)     NATOM
      WRITE(6,2020) NATOM
 2020 FORMAT('# NUMBER OF ATOMS PER CELL=',I4)
      DO 10 I = 0, NATOM-1
         READ(5,*)     XX(I), YY(I)
         WRITE(6,2030) XX(I), YY(I)
 2030 FORMAT('# ',2F6.2)
   10 CONTINUE
C
C     ***** SET ATOM COORDINATES *****
C
      J = 0
      DO 20 IY = 0, NLAT-1
      DO 20 IX = 0, NLAT-1
      DO 20 IN = 0, NATOM - 1
         J = J + 1
         RX(J) = AA*(DFLOAT(IX) + XX(IN))
         RY(J) = AA*(DFLOAT(IY) + YY(IN))
   20 CONTINUE
      JN = J
      WRITE(6,2040) JN
 2040 FORMAT('# JN=',I4)
C
C     ***** SET K-SPACE MESH COORDINATES *****
C
      GG = 2.D0*PI/AA
      DP = 1.D0/DFLOAT(NLAT)
      PMIN = -DFLOAT(NG)
      NP = 2*NG*NLAT
      DO 30 IY = 0, NP
      DO 30 IX = 0, NP
         PX(IX,IY) = PMIN + DP*DFLOAT(IX)
         PY(IX,IY) = PMIN + DP*DFLOAT(IY)
   30 CONTINUE
C
C     ***** CALCULATE STRUCTURE FACTOR *****
C
      DNN = 1.D0/DFLOAT(JN)
      DO 40 JY = 0, NP
      DO 40 JX = 0, NP
         PPX = GG*PX(JX, JY)
         PPY = GG*PY(JX, JY)
         WRE = 0.D0
         WIM = 0.D0
         DO 50 J = 1, JN
            RRX = RX(J)
            RRY = RY(J)
            WRE = WRE + DCOS(PPX*RRX+PPY*RRY)
            WIM = WIM + DSIN(PPX*RRX+PPY*RRY)
   50    CONTINUE
         FF(JX, JY) = (WRE*WRE + WIM*WIM)*DNN
   40 CONTINUE
C
C     ***** OUTPUT RESULTS *****
C
      DO 60 JY = 0, NP
      DO 60 JX = 0, NP
         WRITE(6,1010) PX(JX, JY), PY(JX, JY), FF(JX, JY)
   60 CONTINUE
 1010 FORMAT(3E15.6)
      STOP
      END

入力ファイル例：

12  2.D0  2 
1
0.D0  0.D0