← 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_vel.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_vel.F
 
1
      subroutine smd_vel_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_vel,sp_atom
 
11
      character*32 tag,pname
 
12
      logical result
 
13
 
 
14
      pname = "smd_vel_init_system"
 
15
 
 
16
 
 
17
 
 
18
 
 
19
      write(*,*) "in "//pname
 
20
c
 
21
      tag = "atom"
 
22
      call smd_system_get_component(sp_atom,tag,result)
 
23
      if(.not.result)
 
24
     >  call errquit(
 
25
     >       pname//'no component '//tag,0,0)
 
26
 
 
27
 
 
28
      tag = "velocity"
 
29
      call smd_system_get_component(sp_vel,tag,result)
 
30
      if(.not.result)
 
31
     >  call errquit(
 
32
     >       pname//'no component '//tag,0,0)
 
33
 
 
34
      call smd_vel_init(sp_vel)
 
35
c
 
36
      return
 
37
      end
 
38
 
 
39
      subroutine smd_vel_init(sp_vel)
 
40
      implicit none
 
41
#include "errquit.fh"
 
42
#include "inp.fh"
 
43
#include "mafdecls.fh"
 
44
#include "rtdb.fh"
 
45
#include "util.fh"
 
46
#include "global.fh"
 
47
c     
 
48
      character*(*) sp_vel
 
49
c
 
50
      character*32 pname
 
51
      integer na
 
52
      character*255 filename
 
53
      logical result
 
54
c
 
55
      pname = "smd_vel_init"
 
56
c
 
57
c      write(*,*) "in "//pname
 
58
c
 
59
c     get total number of atoms 
 
60
c     -------------------------
 
61
      call smd_atom_ntot(na)
 
62
      if(na.le.0)
 
63
     >  call errquit(
 
64
     >       pname//'no atoms ',0, RTDB_ERR)
 
65
c
 
66
c     create vel data structures
 
67
c     ---------------------------
 
68
      call smd_namespace_create(sp_vel)
 
69
      call smd_data_create(sp_vel,"vel",3*na,MT_DBL)
 
70
c
 
71
      call smd_velocfile_input(filename,result)
 
72
c      call smd_rtdb_get_veloc_input(filename,result)
 
73
      if((.not.result).or.(filename.eq."random")) then
 
74
        call smd_vel_random_guess()
 
75
      else
 
76
        call smd_vel_read(filename)
 
77
      end if
 
78
c
 
79
      return
 
80
      end
 
81
c
 
82
      subroutine smd_vel_random_guess()
 
83
      implicit none
 
84
#include "errquit.fh"
 
85
#include "inp.fh"
 
86
#include "mafdecls.fh"
 
87
#include "rtdb.fh"
 
88
#include "util.fh"
 
89
c     
 
90
      character*32 pname
 
91
      character*80 tag
 
92
      integer nt,na
 
93
      integer i_v,i_it
 
94
      integer i_m
 
95
      logical result
 
96
      double precision targetke,ke
 
97
c
 
98
      pname = "smd_type_init"
 
99
c
 
100
      write(*,*) "in "//pname
 
101
 
 
102
c
 
103
c     get velocity array
 
104
c     ------------------
 
105
      tag = "vel"
 
106
      call smd_get_ind_dim(tag,i_v,na,result)
 
107
      if(.not. result) 
 
108
     >  call errquit(
 
109
     >       pname//'error getting index for'//tag,0, 0)
 
110
      na = na/3
 
111
 
 
112
c
 
113
c     get mass array
 
114
c     ------------------
 
115
      tag = "param:mass"
 
116
      call smd_get_ind(tag,i_m,result)
 
117
      if(.not. result) 
 
118
     >  call errquit(
 
119
     >       pname//'error getting index for'//tag,0, 0)
 
120
c
 
121
c     get type array
 
122
c     --------------
 
123
      tag = "type:id"
 
124
      call smd_get_ind(tag,i_it,result)
 
125
      if(.not. result) 
 
126
     >  call errquit(
 
127
     >       pname//'error getting index for'//tag,0, 0)
 
128
 
 
129
      call smd_temper_get_ke_target(targetke)
 
130
      if(targetke.eq.0.0d0) 
 
131
     >   call smd_temper_ke_compute(298.15,targetke)
 
132
      call smd_vel_random0(na,
 
133
     >                     targetke,
 
134
     >                     dbl_mb(i_v),
 
135
     >                     dbl_mb(i_m),
 
136
     >                     int_mb(i_it))
 
137
 
 
138
      call smd_vel_ke_compute(ke)
 
139
c      write(*,*) "current ke compared",targetke,ke
 
140
      return
 
141
      end
 
142
 
 
143
      subroutine smd_vel_write(filename)
 
144
      implicit none
 
145
#include "errquit.fh"
 
146
#include "inp.fh"
 
147
#include "mafdecls.fh"
 
148
#include "rtdb.fh"
 
149
#include "util.fh"
 
150
      character*(*) filename
 
151
c
 
152
      integer un
 
153
      character*72 sp_vel
 
154
      character*72 sp_param
 
155
      character*72 sp_type
 
156
      character*32 pname
 
157
      character*80 tag
 
158
      integer nt,na
 
159
      integer i_v,i_it
 
160
      integer i_m
 
161
      logical result
 
162
      double precision targetke,ke
 
163
c
 
164
      pname = "smd_type_init"
 
165
c
 
166
c      write(*,*) "in "//pname
 
167
c
 
168
      if(.not.util_get_io_unit(un))
 
169
     >   call errquit("cannot get file number",0,0)
 
170
c
 
171
      open(unit=un,status="old",form="formatted",file=filename)
 
172
c
 
173
c     get velocity array
 
174
c     ------------------
 
175
      tag = "vel"
 
176
      call smd_get_ind_dim(sp_vel,tag,i_v,na,result)
 
177
      if(.not. result) 
 
178
     >  call errquit(
 
179
     >       pname//'error getting index for'//tag,0, 0)
 
180
      na = na/3
 
181
 
 
182
      call smd_vel_write0(un,
 
183
     >                   na,
 
184
     >                   dbl_mb(i_v))
 
185
 
 
186
      close(un)
 
187
      return
 
188
      end
 
189
 
 
190
      subroutine smd_vel_read(filename)
 
191
      implicit none
 
192
#include "errquit.fh"
 
193
#include "inp.fh"
 
194
#include "mafdecls.fh"
 
195
#include "rtdb.fh"
 
196
#include "util.fh"
 
197
      character*(*) filename
 
198
c
 
199
      integer un
 
200
      character*72 sp_vel
 
201
      character*72 sp_param
 
202
      character*72 sp_type
 
203
      character*32 pname
 
204
      character*80 tag
 
205
      integer nt,na
 
206
      integer i_v,i_it,n
 
207
      integer i_m
 
208
      logical result
 
209
c
 
210
      pname = "smd_vel_read"
 
211
c
 
212
c      write(*,*) "in "//pname
 
213
c
 
214
      if(.not.util_get_io_unit(un))
 
215
     >   call errquit("cannot get file number",0,0)
 
216
c
 
217
      open(unit=un,status="old",form="formatted",file=filename)
 
218
 
 
219
      call smd_get_ind_dim("vel",i_v,na,result)
 
220
      if(.not. result)
 
221
     >  call errquit(
 
222
     >       pname//'error getting index for'//tag,0, 0)
 
223
      na=na/3
 
224
      call smd_vel_read0(un,
 
225
     >                   na,
 
226
     >                   dbl_mb(i_v))
 
227
 
 
228
      close(un)
 
229
      return
 
230
      end
 
231
 
 
232
      subroutine smd_vel_ke_compute(targetke)
 
233
      implicit none
 
234
#include "errquit.fh"
 
235
#include "inp.fh"
 
236
#include "mafdecls.fh"
 
237
#include "rtdb.fh"
 
238
#include "util.fh"
 
239
c     
 
240
      character*32 pname
 
241
      character*80 tag
 
242
      integer nt,na
 
243
      integer i_v,i_it
 
244
      integer i_m
 
245
      logical result
 
246
      double precision targetke
 
247
c
 
248
      pname = "smd_type_init"
 
249
c
 
250
c      write(*,*) "in "//pname
 
251
c
 
252
 
 
253
c
 
254
c     get velocity array
 
255
c     ------------------
 
256
      tag = "vel"
 
257
      call smd_get_ind_dim(tag,i_v,na,result)
 
258
      if(.not. result) 
 
259
     >  call errquit(
 
260
     >       pname//'error getting index for'//tag,0, 0)
 
261
      na = na/3
 
262
 
 
263
c
 
264
c     get mass array
 
265
c     ------------------
 
266
      tag = "param:mass"
 
267
      call smd_get_ind(tag,i_m,result)
 
268
      if(.not. result) 
 
269
     >  call errquit(
 
270
     >       pname//'error getting index for'//tag,0, 0)
 
271
c
 
272
c     get type array
 
273
c     --------------
 
274
      tag = "type:id"
 
275
      call smd_get_ind(tag,i_it,result)
 
276
      if(.not. result) 
 
277
     >  call errquit(
 
278
     >       pname//'error getting index for'//tag,0, 0)
 
279
 
 
280
      call smd_vel_ke0(na,
 
281
     >                     targetke,
 
282
     >                     dbl_mb(i_v),
 
283
     >                     dbl_mb(i_m),
 
284
     >                     int_mb(i_it))
 
285
 
 
286
      call smd_energy_set_component("kinetic",targetke)
 
287
      write(*,*) "current ke 1",targetke
 
288
      return
 
289
      end
 
290
 
 
291
      SUBROUTINE smd_vel_random0(natms,targetke,vvv,typmass,atmtype)
 
292
 
 
293
      implicit none
 
294
 
 
295
      integer natms
 
296
      double precision targetke
 
297
      double precision vvv(natms,3)
 
298
      double precision typmass(*)
 
299
      integer atmtype(*)
 
300
c
 
301
      integer i,iatm,iseed
 
302
 
 
303
      double precision   x
 
304
      double precision  commass,comxvv,comyvv,comzvv
 
305
      double precision  instanke,xscale
 
306
 
 
307
      iseed=620419483
 
308
      comxvv=0.d0
 
309
      comyvv=0.d0
 
310
      comzvv=0.d0
 
311
      commass = 0.0d0
 
312
 
 
313
      do i=1,natms
 
314
 
 
315
       iatm=atmtype(i)
 
316
       call tool_randm(iseed,x)
 
317
       vvv(i,1)=(x-0.5)/sqrt(typmass(iatm))
 
318
       call tool_randm(iseed,x)
 
319
       vvv(i,2)=(x-0.5)/sqrt(typmass(iatm))
 
320
       call tool_randm(iseed,x)
 
321
       vvv(i,3)=(x-0.5)/sqrt(typmass(iatm))
 
322
 
 
323
       comxvv=comxvv+vvv(i,1)*typmass(iatm)
 
324
       comyvv=comyvv+vvv(i,2)*typmass(iatm)
 
325
       comzvv=comzvv+vvv(i,3)*typmass(iatm)
 
326
       commass=commass+typmass(iatm)
 
327
 
 
328
      enddo
 
329
 
 
330
      comxvv=comxvv/commass
 
331
      comyvv=comyvv/commass
 
332
      comzvv=comzvv/commass
 
333
 
 
334
      do i=1,natms
 
335
 
 
336
       vvv(i,1)=vvv(i,1)-comxvv
 
337
       vvv(i,2)=vvv(i,2)-comyvv
 
338
       vvv(i,3)=vvv(i,3)-comzvv
 
339
 
 
340
      enddo
 
341
 
 
342
      instanke=0.d0
 
343
 
 
344
      do i=1,natms
 
345
 
 
346
       iatm=atmtype(i)
 
347
       instanke=instanke
 
348
     $         +typmass(iatm)*((vvv(i,1)**2+vvv(i,2)**2+vvv(i,3)**2))
 
349
 
 
350
      enddo
 
351
 
 
352
      instanke=0.5*instanke
 
353
 
 
354
      xscale=sqrt(targetke/instanke)
 
355
 
 
356
      do i=1,natms
 
357
 
 
358
       vvv(i,1)=xscale*vvv(i,1)
 
359
       vvv(i,2)=xscale*vvv(i,2)
 
360
       vvv(i,3)=xscale*vvv(i,3)
 
361
 
 
362
      enddo
 
363
 
 
364
      instanke=0.d0
 
365
 
 
366
      do i=1,natms
 
367
 
 
368
       iatm=atmtype(i)
 
369
       instanke=instanke
 
370
     $         +typmass(iatm)*((vvv(i,1)**2+vvv(i,2)**2+vvv(i,3)**2))
 
371
 
 
372
      enddo
 
373
 
 
374
      instanke=0.5*instanke
 
375
 
 
376
c      write(*,*) "current ke comp1",targetke,instanke
 
377
 
 
378
      return
 
379
 
 
380
      END
 
381
 
 
382
      SUBROUTINE smd_vel_read0(un,natms,vvv)
 
383
 
 
384
      implicit none
 
385
 
 
386
      integer un,natms
 
387
      double precision vvv(natms,3)
 
388
c
 
389
      integer i
 
390
 
 
391
      do i=1,natms
 
392
 
 
393
       read(un,*) vvv(i,1),vvv(i,2),vvv(i,3)
 
394
 
 
395
      enddo
 
396
 
 
397
      return
 
398
 
 
399
      END
 
400
 
 
401
      SUBROUTINE smd_vel_write0(un,natms,vvv)
 
402
 
 
403
      implicit none
 
404
 
 
405
      integer un,natms
 
406
      double precision vvv(natms,3)
 
407
c
 
408
      integer i
 
409
 
 
410
      do i=1,natms
 
411
 
 
412
       write(un,*) vvv(i,1),vvv(i,2),vvv(i,3)
 
413
 
 
414
      enddo
 
415
 
 
416
      return
 
417
 
 
418
      END
 
419
 
 
420
      SUBROUTINE smd_vel_ke0(natms,ke,vvv,typmass,atmtype)
 
421
 
 
422
      implicit none
 
423
 
 
424
      integer natms
 
425
      double precision targetke
 
426
      double precision vvv(natms,3)
 
427
      double precision typmass(*)
 
428
      integer atmtype(*)
 
429
c
 
430
      integer i,iatm
 
431
 
 
432
      double precision ke
 
433
c
 
434
      ke=0.d0
 
435
 
 
436
      do i=1,natms
 
437
 
 
438
       iatm=atmtype(i)
 
439
       ke=ke
 
440
     $          +typmass(iatm)*((vvv(i,1)**2+vvv(i,2)**2+vvv(i,3)**2))
 
441
 
 
442
      enddo
 
443
 
 
444
      ke=0.5*ke
 
445
 
 
446
      write(*,*) "kinetic energy",ke
 
447
      return
 
448
 
 
449
      END
 
450
 
 
451
      subroutine smd_velocfile_input(filename,result)
 
452
      implicit none
 
453
#include "rtdb.fh"
 
454
#include "smd_rtdb_data.fh"
 
455
#include "mafdecls.fh"
 
456
#include "global.fh"
 
457
 
 
458
c
 
459
      character*(*) filename
 
460
      logical result
 
461
c
 
462
      character*30 pname
 
463
 
 
464
      pname = "smd_velocfile_input"
 
465
 
 
466
      result = .true.
 
467
      call smd_rtdb_get_string("smd:veloc:input",1,
 
468
     >                           filename,result)
 
469
 
 
470
      end