← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/tools/ga-5-2/global/examples/md_cluster/rdpar.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-2/global/examples/md_cluster/rdpar.F
 
1
#if HAVE_CONFIG_H
 
2
#   include "config.fh"
 
3
#endif
 
4
c
 
5
c                                   NOTICE
 
6
c
 
7
c   This software is being made available for internal testing and
 
8
c   evaluation purposes only. This software is a pre-release test version
 
9
c   which has not yet been authenticated and cleared for publication. Adherence
 
10
c   to this notice may be necessary for the author, Battelle Memorial
 
11
c   Institute, to successfully assert copyright in and commercialize this
 
12
c   software. This software is not intended for duplication or distribution
 
13
c   to third parties without the permission of the Manager of Software
 
14
c   Products at Pacific Northwest Laboratory, Richland, Washington,  99352.
 
15
c
 
16
      subroutine rdpar
 
17
#include "common.fh"
 
18
      double precision x,ran1
 
19
      integer i,j
 
20
      integer me
 
21
c
 
22
c    This subroutine reads in the molecular simulation parameters
 
23
c
 
24
      me = ga_nodeid()
 
25
      if (me.eq.0) open (unit=5,file='md_lj.in',status='old')
 
26
c
 
27
c     Line 1
 
28
c
 
29
c       nstep: total number of steps in simulation
 
30
c
 
31
c       tau: time step increment
 
32
c
 
33
c       nsc: number of special control steps
 
34
c
 
35
c       iseed: random number seed (this should be a NEGATIVE
 
36
c              integer)
 
37
c
 
38
      nstep = 0
 
39
      tau = 0.0d00
 
40
      nsc = 0
 
41
      iseed = 0
 
42
      if (me.eq.0) then
 
43
        read(5,*) nstep,tau,nsc,iseed
 
44
      endif
 
45
      call ga_igop(2,nstep,1,'+')
 
46
      call ga_dgop(3,tau,1,'+')
 
47
      call ga_igop(4,nsc,1,'+')
 
48
      call ga_igop(5,iseed,1,'+')
 
49
c
 
50
c     Line 2
 
51
c
 
52
c       rcut: cutoff for use in neighbor list calculation
 
53
c
 
54
c       ilist: number of steps taken before updating
 
55
c              neighbor list
 
56
c
 
57
c       icut: cutoff parameter for potential
 
58
c
 
59
c           (0) only potential set to zero at cutoff
 
60
c
 
61
c           (1) potential and forces set to zero at cutoff
 
62
c
 
63
      rcut = 0.0d00
 
64
      ilist = 0
 
65
      icut = 0
 
66
      if (me.eq.0) then
 
67
        read(5,*) rcut,ilist,icut
 
68
      endif
 
69
      call ga_dgop(6,rcut,1,'+')
 
70
      call ga_igop(7,ilist,1,'+')
 
71
      call ga_igop(8,icut,1,'+')
 
72
c
 
73
c     Line 3
 
74
c
 
75
c       dflalg: default algorithm
 
76
c
 
77
c            1: constant energy 3 point Gear algorithm
 
78
c            2: constant pressure algorithm
 
79
c            3: constant temperature algorithm (velocity scaling)
 
80
c            4: constant temperature algorithm (Nose' dynmamics)
 
81
c            5: constant pressure and temperature algorithm
 
82
c            6: Boltzmann velocity step
 
83
c            7: adjust volume to specified target volume
 
84
c               (available only as a special step)
 
85
c            8: adjust kinetic energy based on average temperature
 
86
c               since last call to this step so that temperature
 
87
c               matches target temperature.
 
88
c               (available only as a special step)
 
89
c            9: constant pressure and temperature algorithm with
 
90
c               independent box adjustments
 
91
c            10: constant pressure and temperature algorithm with
 
92
c               independent box adjustments in only the x and y
 
93
c               dimensions
 
94
c
 
95
c       dftmp: default temperature
 
96
c
 
97
c       dfprs: default pressure
 
98
c
 
99
c       dftm: default mass for temperature
 
100
c
 
101
c       dfpm: default mass for pressure
 
102
c
 
103
      dflalg = 0
 
104
      dftmp = 0.0d00
 
105
      dfprs = 0.0d00
 
106
      dftm = 0.0d00
 
107
      dfpm = 0.0d00
 
108
      if (me.eq.0) then
 
109
        read(5,*) dflalg,dftmp,dfprs,dftm,dfpm
 
110
      endif
 
111
      call ga_igop(9,dflalg,1,'+')
 
112
      call ga_dgop(1,dftmp,1,'+')
 
113
      call ga_dgop(2,dfprs,1,'+')
 
114
      call ga_dgop(3,dftm,1,'+')
 
115
      call ga_dgop(4,dfpm,1,'+')
 
116
c
 
117
c     Line 4.1 - 4.nsc
 
118
c
 
119
c     special step control instructions
 
120
c
 
121
c       isc(i,1): begining step for special instruction
 
122
c
 
123
c       isc(i,2): final step for special instruction
 
124
c
 
125
c       isc(i,3): increment for special instruction
 
126
c
 
127
c       isc(i,4): stepping algorithm for special instruction
 
128
c                 (see documentation on dflalg)
 
129
c
 
130
c       rsc(i,1): temperature (in K) if temperature scaling or
 
131
c                 Boltzmann algorithm is specified
 
132
c
 
133
c       rsc(i,2): pressure if constant pressure
 
134
c                 algorithms are specified. If volume adjustment
 
135
c                 algorithm is specified then this is the target
 
136
c                 volume in A**3.
 
137
c
 
138
c       rsc(i,3): mass (in g A**2 /mole) for Nose temperature
 
139
c                 algorithm
 
140
c
 
141
c       rsc(i,4): mass (in g/(mole A**4)) for constant pressure
 
142
c                 algorithm
 
143
c
 
144
      do 100 i=1,nsc
 
145
        do j=1,4 
 
146
          isc(i,j) = 0.0d00
 
147
          rsc(i,j) = 0.0d00
 
148
        end do
 
149
        if (me.eq.0) then
 
150
          read(5,*) (isc(i,j),j=1,4),(rsc(i,j),j=1,4)
 
151
        endif
 
152
        do j=1,4 
 
153
          call ga_igop(4,isc(i,j),1,'+')
 
154
          call ga_dgop(5,rsc(i,j),1,'+')
 
155
        end do
 
156
  100 continue
 
157
c
 
158
c     Line 5
 
159
c
 
160
c       istart: startup format
 
161
c               1: start from single configuration
 
162
c               2: start from configuration plus velocity
 
163
c
 
164
c       istop: write-out format
 
165
c               1: single configuration
 
166
c               2: configuration plus velocity
 
167
c
 
168
      istart = 0
 
169
      istop = 0
 
170
      if (me.eq.0) then
 
171
        read(5,*) istart,istop
 
172
      endif
 
173
      call ga_igop(6,istart,1,'+')
 
174
      call ga_igop(7,istop,1,'+')
 
175
c
 
176
c     Line 6
 
177
c
 
178
c     istat: frequency to print out simulation information 
 
179
c
 
180
      istat = 0
 
181
      if (me.eq.0) then
 
182
        read(5,*) istat
 
183
      endif
 
184
      call ga_igop(8,istat,1,'+')
 
185
c
 
186
c     Line 7
 
187
c
 
188
c     equil_1: number of steps in first equilibration regime
 
189
c
 
190
c     equil_2: number of steps in second equilibration regime
 
191
c
 
192
c     equil_3: number of steps in third equilibration regime
 
193
c
 
194
      equil_1 = 0
 
195
      equil_2 = 0
 
196
      equil_3 = 0
 
197
      if (me.eq.0) then
 
198
        read(5,*) equil_1, equil_2, equil_3
 
199
      endif
 
200
      call ga_igop(9,equil_1,1,'+')
 
201
      call ga_igop(1,equil_2,1,'+')
 
202
      call ga_igop(2,equil_3,1,'+')
 
203
      window_1 = int(0.2d00*dble(equil_3-equil_2))
 
204
      window_1 = equil_2 + window_1
 
205
      window_2 = int(0.4d00*dble(equil_3-equil_2))
 
206
      window_2 = window_1 + window_2
 
207
c
 
208
c     Line 7
 
209
c
 
210
c     cl_prssr: pressure that is applied to cluster
 
211
c
 
212
c     mc_tmprtr: temperature used in Monte Carlo steps
 
213
c
 
214
      cl_prssr = 0.0d00
 
215
      mc_tmprtr = 0.0d00
 
216
      if (me.eq.0) then
 
217
        read(5,*) cl_prssr, mc_tmprtr
 
218
      endif
 
219
      call ga_dgop(3,cl_prssr,1,'+')
 
220
      call ga_dgop(4,mc_tmprtr,1,'+')
 
221
c
 
222
c     Line 8
 
223
c
 
224
c     mcfreq: frequency to apply Monte Carlo step on radius
 
225
c
 
226
c     mcbins: number of bins to accumulate radius in
 
227
c
 
228
      mcfreq = 0
 
229
      mcbins = 0
 
230
      if (me.eq.0) then
 
231
        read(5,*) mcfreq, mcbins
 
232
      endif
 
233
      call ga_igop(7,mcfreq,1,'+')
 
234
      call ga_igop(8,mcbins,1,'+')
 
235
c
 
236
      if (me.eq.0) close(5)
 
237
c
 
238
c    Write out simulation information to output file
 
239
c
 
240
      if (me.eq.0.and.l_stdio) then
 
241
        write(6,1300) nstep
 
242
        write(6,1400) 1000.0 * tau
 
243
        write(6,4325) rcut
 
244
        write(6,4350) ilist
 
245
        if (icut.eq.0) then
 
246
          write(6,4400)
 
247
        else
 
248
          write(6,4420)
 
249
        endif
 
250
        write(6,1500) dflalg
 
251
        write(6,1600)
 
252
        write(6,1700) dftmp
 
253
        write(6,1800) dfprs
 
254
        write(6,1900) dftm
 
255
        write(6,2000) dfpm
 
256
        write(6,2100)
 
257
        do 300 i = 1, nsc
 
258
          if (i.eq.1) then
 
259
             write(6,2200)
 
260
          else
 
261
             if (isc(i-1,4).eq.7) write(6,2200)
 
262
          endif
 
263
          if (isc(i,4).ne.7) then
 
264
            write(6,2300) (isc(i,j),j=1,4),(rsc(i,j),j=1,4)
 
265
          else
 
266
            write(6,2400) 
 
267
            write(6,2500) (isc(i,j),j=1,4)
 
268
            write(6,2600) rsc(i,1)
 
269
            write(6,2700) rsc(i,2)
 
270
          endif
 
271
  300   continue
 
272
        write(6,3400) istat
 
273
        write(6,3800) mcfreq
 
274
        write(6,3500) equil_1
 
275
        write(6,3600) equil_2
 
276
        write(6,3700) equil_3
 
277
        if (istart.eq.1) then
 
278
          write(6,2800)
 
279
        else
 
280
          write(6,2900)
 
281
        endif
 
282
        if (istop.eq.1) then
 
283
          write(6,3000)
 
284
        else
 
285
          write(6,3100)
 
286
        endif
 
287
        write(6,3200) iseed
 
288
        write(6,4500) ga_pgroup_nodeid(ga_pgroup_get_world())
 
289
      endif
 
290
      iseed = iseed - 10*me
 
291
      x = ran1(iseed)
 
292
      r_cluster = 0.0d00
 
293
      mc_cnt = 0
 
294
c
 
295
      return
 
296
 1300 format('Total number of steps in simulation                :',i9)
 
297
 1400 format('Time step interval                                 :'
 
298
     +       ,f16.6)
 
299
 1500 format('Default algorithm                                  :',i9)
 
300
 1600 format('     Default algorithm parameters ')
 
301
 1700 format('        Default temperature                        :'
 
302
     +       ,f16.6)
 
303
 1800 format('        Default pressure                           :'
 
304
     +       ,f16.6)
 
305
 1900 format('        Default temperature mass                   :'
 
306
     +       ,f16.6)
 
307
 2000 format('        Default pressure mass                      :'
 
308
     +       ,f16.6)
 
309
 2100 format('Special step instructions:')
 
310
 2200 format('   begin     end  incrmt algrthm    T       P      ',
 
311
     +'T mass  P mass  ')
 
312
 2300 format(4i8,2f8.3,1pe9.2,1pe9.2)
 
313
 2400 format('   begin     end  incrmt algrthm')
 
314
 2500 format(4i8)
 
315
 2600 format('     Temperature for volume adjustment             :'
 
316
     +       ,f16.6)
 
317
 2700 format('     Target volume for volume adjustment           :'
 
318
     +       ,f16.6)
 
319
 2800 format('Initial configuration is coordinates only')
 
320
 2900 format('Initial configuration is coordinates and velocities')
 
321
 3000 format('Final configuration is coordinates only')
 
322
 3100 format('Final configuration is coordinates and velocities')
 
323
 3200 format('Random number seed                                 :',i9)
 
324
 3400 format('Frequency to print out simulation information      :',i9)
 
325
 3500 format('Number of steps in first equilibration regime      :',i9)
 
326
 3600 format('Number of steps in second equilibration regime     :',i9)
 
327
 3700 format('Number of steps in third equilibration regime      :',i9)
 
328
 3800 format('Frequency to take Monte Carlo step                 :',i9)
 
329
 4325 format('Cutoff distance                                    :',
 
330
     +       f16.6)
 
331
 4350 format('Neighbor list update frequency                     :',i9)
 
332
 4400 format('Potential set to zero at cutoff distance')
 
333
 4420 format('Potential and forces set to zero at cutoff distance')
 
334
 4500 format('Writing from node                                  :',i9)
 
335
      end
 
336
c
 
337
      subroutine rdcfg
 
338
#include "common.fh"
 
339
c
 
340
      double precision scl,gasdev
 
341
      integer inode,icnt,i,j,k,pnum,ilast,ifirst,me
 
342
      character*32 filename
 
343
c
 
344
c   This subroutine reads in the initial coordinates of all atoms in
 
345
c   the system. First determine which coordinates each processor
 
346
c   should read.
 
347
c
 
348
      me = ga_nodeid()
 
349
      pnum = ga_nnodes()
 
350
c
 
351
c  clean up everything
 
352
c
 
353
      atot = 0
 
354
      ctot = 0
 
355
      antot = 0
 
356
      do k = 1, 8
 
357
        do j = 1, 3
 
358
          do i = 1, MAXAT
 
359
            ra(i,j,k) = 0.0d00
 
360
          end do
 
361
        end do
 
362
      end do
 
363
      do i = 1, MAXAT
 
364
        mass(i) = 0.0d00
 
365
        at(i) = 0
 
366
        aidx(i) = 0
 
367
c
 
368
        xcrd(i) = 0.0d00
 
369
        ycrd(i) = 0.0d00
 
370
        zcrd(i) = 0.0d00
 
371
        xfrc(i) = 0.0d00
 
372
        yfrc(i) = 0.0d00
 
373
        zfrc(i) = 0.0d00
 
374
        xacc(i) = 0.0d00
 
375
        yacc(i) = 0.0d00
 
376
        zacc(i) = 0.0d00
 
377
        mbuf(i) = 0.0d00
 
378
        bat(i) = 0
 
379
        bidx(i) = 0
 
380
      end do
 
381
      btot = 0
 
382
c
 
383
      if (task_id.lt.10) then
 
384
        write(filename,100) task_id
 
385
      else if (task_id.ge.10.and.task_id.lt.100) then
 
386
        write(filename,101) task_id
 
387
      else if (task_id.ge.100.and.task_id.lt.1000) then
 
388
        write(filename,102) task_id
 
389
      else if (task_id.ge.1000.and.task_id.lt.10000) then
 
390
        write(filename,103) task_id
 
391
      endif
 
392
  100 format('md.cfg',i1)
 
393
  101 format('md.cfg',i2)
 
394
  102 format('md.cfg',i3)
 
395
  103 format('md.cfg',i4)
 
396
      atot = 0
 
397
      if (me.eq.0) then
 
398
        open (unit=2,file=filename,status='old')
 
399
        read(2,*) atot
 
400
        close(2)
 
401
      endif
 
402
      call ga_igop(5,atot,1,'+')
 
403
      ilast = nint(dble((me+1)*atot)/dble(pnum))
 
404
      ifirst = nint(dble(me*atot)/dble(pnum))
 
405
      ifirst = ifirst + 1
 
406
c
 
407
c   read in portions of MD configuration on each node
 
408
c
 
409
      do inode = 0, pnum - 1
 
410
        call ga_sync()
 
411
        if (me.eq.inode) then
 
412
          open (unit=2,file=filename,status='old')
 
413
c
 
414
          read(2,*) atot
 
415
          read(2,*) xbox,ybox,zbox
 
416
          xbox2 = xbox/2.0d00
 
417
          ybox2 = ybox/2.0d00
 
418
          zbox2 = zbox/2.0d00
 
419
c
 
420
          cl_lower = 0.0d00
 
421
          if (xbox.lt.ybox.and.xbox.lt.zbox) then
 
422
            cl_upper = xbox2
 
423
          else if (ybox.lt.xbox.and.ybox.lt.zbox) then
 
424
            cl_upper = ybox2
 
425
          else
 
426
            cl_upper = zbox2
 
427
          endif
 
428
          mc_step = (cl_upper-cl_lower)/dble(mcbins)
 
429
c
 
430
          icnt = 0
 
431
          if (istart.eq.1.or.(.not.l_oldcfg)) then
 
432
            do i = 1, ilast
 
433
              if (i.lt.ifirst) then
 
434
                read(2,*)
 
435
              else
 
436
                icnt = icnt + 1
 
437
                read(2,*) at(icnt),(ra(icnt,j,1),j=1,3)
 
438
                if (at(icnt).eq.2) ctot = ctot + 1
 
439
c
 
440
c    generate atomic velocities from a Maxwell-Boltzmann distribution
 
441
c
 
442
                do j = 1, 3
 
443
                  scl = sqrt(dftmp / amass(at(icnt)))
 
444
                  ra(icnt,j,2) = scl * gasdev(0)
 
445
                end do
 
446
                aidx(icnt) = i
 
447
              endif
 
448
            end do
 
449
          else
 
450
            do i = 1, ilast
 
451
              if (i.lt.ifirst) then
 
452
                read(2,*)
 
453
              else
 
454
                icnt = icnt + 1
 
455
                read(2,*) at(icnt),(ra(icnt,j,1),j=1,3),
 
456
     +                    (ra(icnt,j,2),j=1,3)
 
457
                if (at(icnt).eq.2) ctot = ctot + 1
 
458
                aidx(icnt) = i
 
459
              endif
 
460
            end do
 
461
          endif
 
462
          close(2)
 
463
        endif
 
464
      end do
 
465
      antot = icnt
 
466
      call ga_igop(9,ctot,1,'+')
 
467
c
 
468
      if (me.eq.0.and.l_stdio) then
 
469
        write(6,1100) atot
 
470
        write(6,1000) xbox,ybox,zbox
 
471
      endif
 
472
c
 
473
c   initialize absolute coordinates
 
474
c
 
475
      do j = 1, 3
 
476
        do i = 1, antot
 
477
          ra(i,j,6) = ra(i,j,1)
 
478
        end do
 
479
      end do
 
480
c
 
481
      call fixper
 
482
c
 
483
      return
 
484
 1000 format('The initial box size parameters are '/
 
485
     +       '      x dimension = ',f9.4,' A'/
 
486
     +       '      y dimension = ',f9.4,' A'/
 
487
     +       '      z dimension = ',f9.4,' A') 
 
488
 1100 format('The total number of atoms in simulation            :',i9)
 
489
      end
 
490
c
 
491
      subroutine atomin
 
492
#include "common.fh"
 
493
c
 
494
      double precision epsln,sigma
 
495
      integer i,j,ipairs,ip,ktrnc
 
496
      integer me
 
497
      logical ichk(50,50)
 
498
c
 
499
c    This subroutine reads in atomic parameters from a parameter
 
500
c    file. These parameters will be used to construct the
 
501
c    system potential function.
 
502
c
 
503
      me = ga_nodeid()
 
504
      do i = 1, 50
 
505
        do j = 1, 50
 
506
          ichk(i,j) = .false.
 
507
        end do
 
508
      end do
 
509
c
 
510
      if (me.eq.0) then
 
511
        open(unit=2,file='atom.inp',status='old')
 
512
      endif
 
513
c
 
514
      if (me.eq.0) then
 
515
        read (2,*) amass(1)
 
516
        read (2,*) amass(2)
 
517
        read (2,*) amass(3)
 
518
        read (2,*) amass(4)
 
519
        atnum = 4
 
520
      else
 
521
        amass(1) = 0.0d00
 
522
        amass(2) = 0.0d00
 
523
        amass(3) = 0.0d00
 
524
        amass(4) = 0.0d00
 
525
        atnum = 0
 
526
      endif
 
527
      call ga_dgop(1,amass(1),4,'+')
 
528
      call ga_igop(2,atnum,1,'+')
 
529
c
 
530
c   read in atomic parameters
 
531
c
 
532
      if (me.eq.0) then
 
533
        read(2,*) ipairs
 
534
      else
 
535
        ipairs = 0
 
536
      endif
 
537
      call ga_igop(3,ipairs,1,'+')
 
538
c
 
539
      do 100 ip = 1, ipairs
 
540
        i = 0
 
541
        j = 0
 
542
        ktrnc = 0
 
543
        epsln = 0.0d00
 
544
        sigma = 0.0d00
 
545
        if (me.eq.0) then
 
546
          read (2,*) i,j,ktrnc,epsln,sigma
 
547
        endif
 
548
        call ga_igop(1,i,1,'+')
 
549
        call ga_igop(2,j,1,'+')
 
550
        call ga_igop(3,ktrnc,1,'+')
 
551
        call ga_dgop(4,epsln,1,'+')
 
552
        call ga_dgop(5,sigma,1,'+')
 
553
        ichk(i,j) = .true.
 
554
        if (i.ne.j) ichk(j,i) = .true.
 
555
        e12(i,j) = 4.0d00*epsln*sigma**12
 
556
        e6(i,j) = 4.0d00*epsln*sigma**6
 
557
        if (i.ne.j) e12(j,i) = e12(i,j)
 
558
        if (i.ne.j) e6(j,i) = e6(i,j)
 
559
        if (ktrnc.eq.1) then
 
560
          acut(i,j) = exp(log(2.0d00)/6.0d00)*sigma
 
561
          acut(j,i) = exp(log(2.0d00)/6.0d00)*sigma
 
562
        else
 
563
          acut(i,j) = rcut * sigma
 
564
          acut(j,i) = rcut * sigma
 
565
        endif
 
566
        acut2(i,j) = acut(i,j)**2
 
567
        acut2(j,i) = acut(j,i)**2
 
568
  100 continue
 
569
c
 
570
      if (me.eq.0) close(2)
 
571
      return
 
572
      end