← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/smd/graveyard/smd-9-10-08/smd_vlist.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/smd/graveyard/smd-9-10-08/smd_vlist.F
 
1
      subroutine smd_vlist_init_system()
 
2
      implicit none
 
3
#include "errquit.fh"
 
4
#include "inp.fh"
 
5
#include "mafdecls.fh"
 
6
#include "rtdb.fh"
 
7
#include "util.fh"
 
8
#include "global.fh"
 
9
c     
 
10
      character*32 sp_vlist,sp_exlist,sp_coords
 
11
      character*32 tag,pname
 
12
      logical result
 
13
      logical oexlist
 
14
 
 
15
      pname = "smd_vlist_init_system"
 
16
c
 
17
      tag = "coordinates"
 
18
      call smd_system_get_component(sp_coords,tag,result)
 
19
      if(.not.result)
 
20
     >  call errquit(
 
21
     >       pname//'no component '//tag,0,0)
 
22
 
 
23
      oexlist = .true.
 
24
      tag = "excl_list"
 
25
      call smd_system_get_component(sp_exlist,tag,result)
 
26
      if(.not.result) oexlist = .false.
 
27
 
 
28
      tag = "verlet_list"
 
29
      call smd_system_get_component(sp_vlist,tag,result)
 
30
      if(.not.result)
 
31
     >  call errquit(
 
32
     >       pname//'no component '//tag,0,0)
 
33
 
 
34
 
 
35
      call smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coords)
 
36
c
 
37
      return
 
38
      end
 
39
 
 
40
      subroutine smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coord)
 
41
      implicit none
 
42
#include "errquit.fh"
 
43
#include "inp.fh"
 
44
#include "mafdecls.fh"
 
45
#include "rtdb.fh"
 
46
#include "util.fh"
 
47
#include "global.fh"
 
48
c     
 
49
      character*(*) sp_vlist
 
50
      character*(*) sp_exlist
 
51
      character*(*) sp_coord
 
52
      logical oexlist
 
53
c
 
54
      character*32 pname
 
55
      character*80 tag
 
56
      integer np,nl
 
57
      integer i_p,i_l,i_c
 
58
      integer h_l
 
59
      integer i_list,i_clist
 
60
      integer i
 
61
      integer i_xp,i_xl
 
62
      integer i_cl,h_cl
 
63
      integer i_ct,h_ct
 
64
      double precision rc2
 
65
      integer nlb
 
66
      logical result
 
67
c
 
68
      pname = "smd_vlist_init"
 
69
c
 
70
c      write(*,*) "in "//pname
 
71
c
 
72
c     get coordinate information
 
73
c     --------------------------
 
74
      tag = "coords"
 
75
      call smd_get_ind_dim(tag,i_c,np,result)
 
76
      if(.not. result) 
 
77
     >  call errquit(
 
78
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
79
      np = np/3
 
80
c
 
81
c     get excluded list information
 
82
c     -----------------------------
 
83
      if(oexlist) then
 
84
      tag = "exlist:pointer"
 
85
      call smd_get_ind(tag,i_xp,result)
 
86
      if(.not. result) 
 
87
     >  call errquit(
 
88
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
89
      tag = "exlist:list"
 
90
      call smd_get_ind(tag,i_xl,result)
 
91
      if(.not. result) 
 
92
     >  call errquit(
 
93
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
94
 
 
95
      end if
 
96
c
 
97
c     gestimate the size of pair list
 
98
c     ------------------------------
 
99
      nl =  min( 7*np*200, ma_inquire_avail(MT_INT))
 
100
      nl = nl/7
 
101
c
 
102
c     create pointer memory
 
103
c     ---------------------
 
104
      call smd_namespace_create(sp_vlist)
 
105
      tag = "vlist:pointer"
 
106
      call smd_data_create_get(sp_vlist,"vlist:pointer",np,MT_INT,i_p)
 
107
c
 
108
c    create temporary scratch array for list since
 
109
c    we do not know the size yet
 
110
c    ---------------------------------------------
 
111
      if(.not.ma_push_get(mt_int,nl,'tmp l',h_l,i_l))
 
112
     + call errquit(pname//'Failed to allocate memory for tmp l',
 
113
     + nl, MA_ERR)
 
114
 
 
115
      if(.not.ma_push_get(mt_dbl,3*nl,'tmp cl',h_cl,i_cl))
 
116
     + call errquit(pname//'Failed to allocate memory for tmp l',
 
117
     + nl, MA_ERR)
 
118
 
 
119
      if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct))
 
120
     + call errquit(pname//'Failed to allocate memory for tmp l',
 
121
     + nl, MA_ERR)
 
122
 
 
123
 
 
124
      call smd_cutoff_get_rcut_verlet(rc2)
 
125
      rc2=rc2*rc2
 
126
 
 
127
      if(oexlist) then
 
128
      call smd_vlist_set(np,
 
129
     +                   nl,
 
130
     +                   rc2,
 
131
     +                   dbl_mb(i_c),
 
132
     +                   int_mb(i_xp),
 
133
     +                   int_mb(i_xl),
 
134
     +                   int_mb(i_p),
 
135
     +                   int_mb(i_l),
 
136
     +                   dbl_mb(i_cl),
 
137
     +                   dbl_mb(i_ct),
 
138
     +                   result)
 
139
 
 
140
      else
 
141
 
 
142
      call smd_vlist_set1(np,
 
143
     +                   nl,
 
144
     +                   rc2,
 
145
     +                   dbl_mb(i_c),
 
146
     +                   int_mb(i_p),
 
147
     +                   int_mb(i_l),
 
148
     +                   dbl_mb(i_cl),
 
149
     +                   dbl_mb(i_ct),
 
150
     +                   result)
 
151
 
 
152
 
 
153
      end if
 
154
c
 
155
c     create list memory
 
156
c     nl now contains the actual size
 
157
c     we will buffer it though to allow for possible expansion
 
158
c     --------------------------------------------------------
 
159
      nlb = 500
 
160
      tag = "vlist:list"
 
161
      call smd_data_create_get(sp_vlist,tag,nl+nlb,MT_INT,i_list)
 
162
 
 
163
      tag = "vlist:distances"
 
164
      call smd_data_create_get(sp_vlist,tag,3*(nl+nlb),MT_DBL,i_clist)
 
165
 
 
166
      tag = "vlist:displacement"
 
167
      call smd_data_create(sp_vlist,tag,3*np,MT_DBL)
 
168
 
 
169
      do i=1,3*nl
 
170
       dbl_mb(i_clist+i-1) = dbl_mb(i_cl+i-1)
 
171
      end do 
 
172
       
 
173
      if(.not.ma_pop_stack(h_ct))
 
174
     & call errquit(pname//'Failed to deallocate stack h_l',nl,
 
175
     &       MA_ERR)
 
176
 
 
177
     
 
178
      if(.not.ma_pop_stack(h_cl))
 
179
     & call errquit(pname//'Failed to deallocate stack h_l',nl,
 
180
     &       MA_ERR)
 
181
 
 
182
 
 
183
      if(.not.ma_pop_stack(h_l))
 
184
     & call errquit(pname//'Failed to deallocate stack h_l',nl,
 
185
     &       MA_ERR)
 
186
 
 
187
 
 
188
      return
 
189
      end
 
190
c
 
191
      subroutine smd_vlist_update(olist,ocoord)
 
192
      implicit none
 
193
#include "errquit.fh"
 
194
#include "inp.fh"
 
195
#include "mafdecls.fh"
 
196
#include "rtdb.fh"
 
197
#include "util.fh"
 
198
#include "global.fh"
 
199
c     
 
200
      logical olist
 
201
      logical ocoord
 
202
c
 
203
      character*32 sp_vlist
 
204
      character*32 sp_exlist
 
205
      character*32 sp_coord
 
206
c
 
207
      character*32 pname
 
208
      character*80 tag
 
209
      integer np,nl
 
210
      integer i_p,i_l,i_c
 
211
      integer h_l
 
212
      integer i_list,i_clist
 
213
      integer i
 
214
      integer i_xp,i_xl
 
215
      double precision rc2
 
216
      integer nlb
 
217
      logical result
 
218
      integer i_ct,h_ct
 
219
      logical oexlist 
 
220
c
 
221
      pname = "smd_vlist_update"
 
222
c
 
223
c      write(*,*) "in "//pname
 
224
c
 
225
      if((.not.olist).and.(.not.ocoord)) return
 
226
c
 
227
c     get components
 
228
c     --------------
 
229
      tag = "coordinates"
 
230
      call smd_system_get_component(sp_coord,tag,result)
 
231
      if(.not.result)
 
232
     >  call errquit(
 
233
     >       pname//'no component '//tag,0,0)
 
234
 
 
235
      oexlist = .true.
 
236
      tag = "excl_list"
 
237
      call smd_system_get_component(sp_exlist,tag,result)
 
238
      if(.not.result) oexlist =.false.
 
239
 
 
240
      tag = "verlet_list"
 
241
      call smd_system_get_component(sp_vlist,tag,result)
 
242
      if(.not.result)
 
243
     >  call errquit(
 
244
     >       pname//'no component '//tag,0,0)
 
245
 
 
246
c
 
247
      if(.not.olist) then
 
248
        call smd_vlist_update_coord(sp_vlist,sp_coord)
 
249
        goto 200
 
250
      end if
 
251
c
 
252
c     get coordinate information
 
253
c     --------------------------
 
254
      tag = "coords"
 
255
      call smd_get_ind_dim(tag,i_c,np,result)
 
256
      if(.not. result) 
 
257
     >  call errquit(
 
258
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
259
      np = np/3 
 
260
c
 
261
c     get excluded list information
 
262
c     -----------------------------
 
263
      if(oexlist) then
 
264
      tag = "exlist:pointer"
 
265
      call smd_get_ind(tag,i_xp,result)
 
266
      if(.not. result) 
 
267
     >  call errquit(
 
268
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
269
      tag = "exlist:list"
 
270
      call smd_get_ind(tag,i_xl,result)
 
271
      if(.not. result) 
 
272
     >  call errquit(
 
273
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
274
 
 
275
      end if
 
276
c
 
277
c     get verlet list information
 
278
c     ---------------------------
 
279
      tag = "vlist:pointer"
 
280
      call smd_get_ind(tag,i_p,result)
 
281
      if(.not. result) 
 
282
     >  call errquit(
 
283
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
284
 
 
285
      tag = "vlist:list"
 
286
      call smd_get_ind_dim(tag,i_list,nl,result)
 
287
      if(.not. result) 
 
288
     >  call errquit(
 
289
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
290
 
 
291
      tag = "vlist:distances"
 
292
      call smd_get_ind(tag,i_clist,result)
 
293
      if(.not. result) 
 
294
     >  call errquit(
 
295
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
296
c
 
297
c    create temporary scratch array for list
 
298
c    ---------------------------------------------
 
299
 
 
300
      if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct))
 
301
     + call errquit(pname//'Failed to allocate memory for tmp l',
 
302
     + nl, MA_ERR)
 
303
 
 
304
 
 
305
      call smd_cutoff_get_rcut_verlet(rc2)
 
306
      rc2=rc2*rc2
 
307
 
 
308
      if(oexlist) then
 
309
      call smd_vlist_set(np,
 
310
     +                   nl,
 
311
     +                   rc2,
 
312
     +                   dbl_mb(i_c),
 
313
     +                   int_mb(i_xp),
 
314
     +                   int_mb(i_xl),
 
315
     +                   int_mb(i_p),
 
316
     +                   int_mb(i_list),
 
317
     +                   dbl_mb(i_clist),
 
318
     +                   dbl_mb(i_ct),
 
319
     +                   result)
 
320
 
 
321
      else
 
322
 
 
323
      call smd_vlist_set1(np,
 
324
     +                   nl,
 
325
     +                   rc2,
 
326
     +                   dbl_mb(i_c),
 
327
     +                   int_mb(i_p),
 
328
     +                   int_mb(i_list),
 
329
     +                   dbl_mb(i_clist),
 
330
     +                   dbl_mb(i_ct),
 
331
     +                   result)
 
332
 
 
333
 
 
334
      end if
 
335
       
 
336
      if(.not.ma_pop_stack(h_ct))
 
337
     & call errquit(pname//'Failed to deallocate stack h_l',nl,
 
338
     &       MA_ERR)
 
339
 
 
340
 
 
341
200   continue
 
342
      return
 
343
      end
 
344
c
 
345
      subroutine smd_vlist_update_coord(sp_vlist,sp_coord)
 
346
      implicit none
 
347
#include "errquit.fh"
 
348
#include "inp.fh"
 
349
#include "mafdecls.fh"
 
350
#include "rtdb.fh"
 
351
#include "util.fh"
 
352
#include "global.fh"
 
353
c     
 
354
      character*(*) sp_vlist
 
355
      character*(*) sp_coord
 
356
c
 
357
      character*32 pname
 
358
      character*80 tag
 
359
      integer np,nl
 
360
      integer i_p,i_c
 
361
      integer i_list,i_clist
 
362
      logical result
 
363
      integer i
 
364
c
 
365
      pname = "smd_vlist_init"
 
366
c
 
367
c      write(*,*) "in "//pname
 
368
c
 
369
c     get coordinate information
 
370
c     --------------------------
 
371
      tag = "coords"
 
372
      call smd_get_ind_dim(tag,i_c,np,result)
 
373
      if(.not. result) 
 
374
     >  call errquit(
 
375
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
376
      np = np/3
 
377
c
 
378
c     get list information
 
379
c     --------------------
 
380
      tag = "vlist:pointer"
 
381
      call smd_get_ind(tag,i_p,result)
 
382
      if(.not. result) 
 
383
     >  call errquit(
 
384
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
385
 
 
386
      tag = "vlist:list"
 
387
      call smd_get_ind(tag,i_list,result)
 
388
      if(.not. result) 
 
389
     >  call errquit(
 
390
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
391
 
 
392
 
 
393
      tag = "vlist:distances"
 
394
      call smd_get_ind_dim(tag,i_clist,nl,result)
 
395
      if(.not. result) 
 
396
     >  call errquit(
 
397
     >       pname//'error getting index for'//tag,0, RTDB_ERR)
 
398
      nl = nl/3
 
399
 
 
400
      call smd_vlist_update_coord0(np,nl,
 
401
     +                   dbl_mb(i_c),
 
402
     +                   int_mb(i_p),
 
403
     +                   int_mb(i_list),
 
404
     +                   dbl_mb(i_clist))
 
405
 
 
406
      call smd_lat_rebox(nl,dbl_mb(i_clist))
 
407
 
 
408
      return
 
409
      end
 
410
c
 
411
      subroutine smd_vlist_update_coord0(np,nl,c,point,
 
412
     >                         list,cl)
 
413
c
 
414
c     update coordinates in verlet pair list
 
415
c     point(i) is an index to list() array
 
416
c     that contains all the pairs of atom i
 
417
c     In other words all the atoms paired with atom i
 
418
c     are contained in list(point(i)),...,list(point(i+1)-1)
 
419
c     np     [in]  number of atoms (which is also size of pointer array)
 
420
c     nl     [in]  list size
 
421
c     c      [in]  coordinates
 
422
c     point  [in] verlet pointer
 
423
c     list   [in] verlet list
 
424
c     cl     [out] list of vectors (ri-rj)
 
425
      implicit none
 
426
#include "errquit.fh"
 
427
      integer np,nl
 
428
      double precision c(np,3)
 
429
      integer point(np)
 
430
      integer list(nl)
 
431
      double precision cl(nl,3)
 
432
c
 
433
      integer i,j
 
434
      integer nlist
 
435
      integer jnab,jbeg,jend
 
436
 
 
437
      character*30 pname
 
438
 
 
439
      pname = "smd_vlist_update_coord0"
 
440
      nlist=0
 
441
      do i=1,np-1
 
442
       jbeg=point(i)
 
443
       jend=point(i+1)-1
 
444
 
 
445
       if(jbeg.le.jend)then
 
446
 
 
447
        do jnab=jbeg,jend
 
448
 
 
449
         j=list(jnab)
 
450
 
 
451
         nlist = nlist + 1
 
452
         if(nlist.gt.nl)
 
453
     >       call errquit(
 
454
     >       pname//'out of bounds',0, RTDB_ERR)
 
455
 
 
456
         cl(nlist,1)=c(i,1)-c(j,1)
 
457
         cl(nlist,2)=c(i,2)-c(j,2)
 
458
         cl(nlist,3)=c(i,3)-c(j,3)
 
459
 
 
460
        enddo
 
461
 
 
462
       end if
 
463
      enddo
 
464
 
 
465
 
 
466
200   continue
 
467
      return
 
468
 
 
469
      END
 
470
c
 
471
      subroutine smd_vlist_set(np,nl,vcutsq,c,xp,xl,point,
 
472
     >                         list,cl,ct,result)
 
473
c
 
474
c     constructs verlet pairt list
 
475
c     point(i) is an index to list() array
 
476
c     that contains all the pairs of atom i
 
477
c     In other words all the atoms paired with atom i
 
478
c     are contained in list(point(i)),...,list(point(i+1)-1)
 
479
c     np     [in]  number of atoms (which is also size of pointer array)
 
480
c     nl     [inout] size of verlet list 
 
481
c     c      [in]  coordinates
 
482
c     xp     [in]  excluded pointer 
 
483
c     xl     [in]  excluded list
 
484
c     point  [out] verlet pointer
 
485
c     list   [out] verlet list
 
486
c     cl     [out] list of vectors (ri-rj)
 
487
c     result [out] status of subroutine
 
488
      implicit none
 
489
#include "errquit.fh"
 
490
      integer np
 
491
      integer nl
 
492
      double precision vcutsq
 
493
      double precision c(np,3)
 
494
      integer xp(np)
 
495
      integer xl(*)
 
496
      integer point(np)
 
497
      integer list(nl)
 
498
      double precision cl(nl,3)
 
499
      double precision ct(np,3)
 
500
      logical result
 
501
c
 
502
      integer i,j,k
 
503
      integer nlist
 
504
      integer eatm
 
505
      double precision rij(3),rijsq,cc(1,3)
 
506
 
 
507
      character*30 pname
 
508
 
 
509
      pname = "smd_vlist_set"
 
510
 
 
511
      result = .false.
 
512
      nlist=0
 
513
 
 
514
      do i=1,np-1
 
515
 
 
516
       point(i)=nlist+1
 
517
       if(xp(i).ne.xp(i+1))eatm=xp(i)
 
518
 
 
519
       k = 0
 
520
       do j=i+1,np
 
521
 
 
522
         k = k +1
 
523
         ct(k,1)=c(i,1)-c(j,1)
 
524
         ct(k,2)=c(i,2)-c(j,2)
 
525
         ct(k,3)=c(i,3)-c(j,3)
 
526
 
 
527
       end do
 
528
 
 
529
 
 
530
       call smd_lat_rebox(np,ct)
 
531
 
 
532
 
 
533
       k = 0
 
534
       do j=i+1,np
 
535
 
 
536
        k = k + 1
 
537
        if((xp(i).ne.xp(i+1)).and.(xl(eatm).eq.j))then
 
538
 
 
539
         eatm=min(eatm+1,(xp(i+1)-1))
 
540
 
 
541
        else
 
542
 
 
543
         rij(1)=ct(k,1)
 
544
         rij(2)=ct(k,2)
 
545
         rij(3)=ct(k,3)
 
546
 
 
547
         rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)
 
548
 
 
549
 
 
550
         if(rijsq.lt.vcutsq)then
 
551
 
 
552
          nlist=nlist+1
 
553
 
 
554
          if(nlist.gt.nl)then
 
555
           result = .false.
 
556
           goto 200
 
557
          endif
 
558
 
 
559
          list(nlist)=j
 
560
          cl(nlist,1)=rij(1)
 
561
          cl(nlist,2)=rij(2)
 
562
          cl(nlist,3)=rij(3)
 
563
 
 
564
         endif
 
565
 
 
566
        endif
 
567
 
 
568
       enddo
 
569
 
 
570
      enddo
 
571
 
 
572
      point(np)=nlist+1
 
573
 
 
574
      nl = nlist
 
575
 
 
576
      result = .true.
 
577
 
 
578
200   continue
 
579
      return
 
580
 
 
581
      END
 
582
 
 
583
      subroutine smd_vlist_set1(np,nl,vcutsq,c,point,
 
584
     >                         list,cl,ct,result)
 
585
c
 
586
c     constructs verlet pairt list
 
587
c     point(i) is an index to list() array
 
588
c     that contains all the pairs of atom i
 
589
c     In other words all the atoms paired with atom i
 
590
c     are contained in list(point(i)),...,list(point(i+1)-1)
 
591
c     np     [in]  number of atoms (which is also size of pointer array)
 
592
c     nl     [inout] size of verlet list 
 
593
c     c      [in]  coordinates
 
594
c     xp     [in]  excluded pointer 
 
595
c     xl     [in]  excluded list
 
596
c     point  [out] verlet pointer
 
597
c     list   [out] verlet list
 
598
c     cl     [out] list of vectors (ri-rj)
 
599
c     result [out] status of subroutine
 
600
      implicit none
 
601
#include "errquit.fh"
 
602
      integer np
 
603
      integer nl
 
604
      double precision vcutsq
 
605
      double precision c(np,3)
 
606
      integer point(np)
 
607
      integer list(nl)
 
608
      double precision cl(nl,3)
 
609
      double precision ct(np,3)
 
610
      logical result
 
611
c
 
612
      integer i,j,k
 
613
      integer nlist
 
614
      integer eatm
 
615
      double precision rij(3),rijsq,cc(1,3)
 
616
 
 
617
      character*30 pname
 
618
 
 
619
      pname = "smd_vlist_set"
 
620
 
 
621
      result = .false.
 
622
      nlist=0
 
623
 
 
624
      do i=1,np-1
 
625
 
 
626
       point(i)=nlist+1
 
627
 
 
628
       k = 0
 
629
       do j=i+1,np
 
630
 
 
631
         k = k +1
 
632
         ct(k,1)=c(i,1)-c(j,1)
 
633
         ct(k,2)=c(i,2)-c(j,2)
 
634
         ct(k,3)=c(i,3)-c(j,3)
 
635
 
 
636
       end do
 
637
 
 
638
 
 
639
       call smd_lat_rebox(np,ct)
 
640
 
 
641
 
 
642
       k = 0
 
643
       do j=i+1,np
 
644
 
 
645
        k = k + 1
 
646
 
 
647
         rij(1)=ct(k,1)
 
648
         rij(2)=ct(k,2)
 
649
         rij(3)=ct(k,3)
 
650
 
 
651
         rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)
 
652
 
 
653
 
 
654
         if(rijsq.lt.vcutsq)then
 
655
 
 
656
          nlist=nlist+1
 
657
 
 
658
          if(nlist.gt.nl)then
 
659
           result = .false.
 
660
           goto 200
 
661
          endif
 
662
 
 
663
          list(nlist)=j
 
664
          cl(nlist,1)=rij(1)
 
665
          cl(nlist,2)=rij(2)
 
666
          cl(nlist,3)=rij(3)
 
667
 
 
668
         endif
 
669
 
 
670
       enddo
 
671
 
 
672
      enddo
 
673
 
 
674
      point(np)=nlist+1
 
675
 
 
676
      nl = nlist
 
677
 
 
678
      result = .true.
 
679
 
 
680
200   continue
 
681
      return
 
682
 
 
683
      END
 
684
 
 
685
      SUBROUTINE smd_vlist_test(lupdate)
 
686
      implicit none
 
687
#include "smd_system.fh"
 
688
#include "mafdecls.fh"
 
689
#include "errquit.fh"
 
690
c 
 
691
      logical lupdate
 
692
c
 
693
      integer na
 
694
      integer i_v,i_vd
 
695
      double precision t
 
696
      character*72 sp_vel
 
697
      character*72 sp_vlist
 
698
      logical result
 
699
      character*72 tag
 
700
      character*30 pname
 
701
      double precision vcut,rcut
 
702
c
 
703
      pname = "smd_vlist_test"
 
704
c
 
705
c     get velocity aray
 
706
c     -----------------
 
707
      tag = "velocity"
 
708
      call smd_system_get_component(sp_vel,tag,result)
 
709
      if(.not.result)
 
710
     >  call errquit(
 
711
     >       pname//'no component '//tag,0,0)
 
712
 
 
713
      tag = "vel"
 
714
      call smd_get_ind_dim(tag,i_v,na,result)
 
715
      if(.not. result) 
 
716
     >  call errquit(
 
717
     >       pname//'error getting index for'//tag,0, 0)
 
718
      na = na/3
 
719
c
 
720
c     get time step
 
721
c     ------------
 
722
      if(.not.smd_system_tstep(t)) 
 
723
     >  call errquit(
 
724
     >       pname//'no time step ',0,0)
 
725
 
 
726
c
 
727
c     get verlet displacement array
 
728
c     -----------------------------
 
729
      tag = "vlist:displacement"
 
730
      call smd_get_ind(tag,i_vd,result)
 
731
      if(.not. result) 
 
732
     >  call errquit(
 
733
     >       pname//'error getting index for'//tag,0, 0)
 
734
 
 
735
       call smd_cutoff_get_rcut(rcut)
 
736
       call smd_cutoff_get_rcut_verlet(vcut)
 
737
 
 
738
       call smd_vlist_test0(na,t,rcut,vcut,
 
739
     >     dbl_mb(i_vd),dbl_mb(i_v),lupdate)
 
740
 
 
741
c      write(*,*) "out",pname
 
742
      return
 
743
 
 
744
      END
 
745
 
 
746
      SUBROUTINE smd_vlist_test0(natms,tstep,rcut,vcut,ivv,vvv,lupdate)
 
747
 
 
748
      implicit none
 
749
c
 
750
      integer natms
 
751
      double precision rcut,vcut
 
752
      double precision tstep
 
753
      logical lupdate
 
754
      double precision ivv(natms,3)
 
755
      double precision vvv(natms,3)
 
756
c
 
757
      integer i,exceed
 
758
 
 
759
      double precision  tstepsq
 
760
      double precision  dispmax,dispxsq,dispysq,dispzsq,disprsq 
 
761
 
 
762
      logical lnew
 
763
 
 
764
      data lnew/.true./
 
765
 
 
766
      save lnew
 
767
 
 
768
      tstepsq=tstep**2
 
769
 
 
770
      if(lnew)then
 
771
 
 
772
       lupdate=.true.
 
773
       lnew=.false.
 
774
 
 
775
       do i=1,natms
 
776
 
 
777
        ivv(i,1)=0.0
 
778
        ivv(i,2)=0.0
 
779
        ivv(i,3)=0.0
 
780
 
 
781
       enddo
 
782
 
 
783
      else
 
784
 
 
785
       lupdate=.false.
 
786
 
 
787
       dispmax=((vcut-rcut)/2.0)**2
 
788
 
 
789
       do i=1,natms
 
790
 
 
791
        ivv(i,1)=ivv(i,1)+vvv(i,1)
 
792
        ivv(i,2)=ivv(i,2)+vvv(i,2)
 
793
        ivv(i,3)=ivv(i,3)+vvv(i,3)
 
794
 
 
795
       enddo
 
796
 
 
797
       exceed=0
 
798
 
 
799
       do i=1,natms
 
800
 
 
801
        dispxsq=ivv(i,1)**2
 
802
        dispysq=ivv(i,2)**2
 
803
        dispzsq=ivv(i,3)**2
 
804
        disprsq=tstepsq*(dispxsq+dispysq+dispzsq)
 
805
        if(disprsq.gt.dispmax) then
 
806
          exceed=exceed+1
 
807
          write(*,*) "verlet update disp",disprsq,dispmax
 
808
        end if
 
809
        if(exceed.ge.2)lupdate=.true.
 
810
       enddo
 
811
 
 
812
       if(lupdate)then
 
813
 
 
814
        do i=1,natms
 
815
 
 
816
         ivv(i,1)=0
 
817
         ivv(i,2)=0
 
818
         ivv(i,3)=0
 
819
 
 
820
        enddo
 
821
 
 
822
       endif
 
823
 
 
824
      endif
 
825
 
 
826
      return
 
827
 
 
828
      END