← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/optim/neb/bead_list.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/optim/neb/bead_list.F
1
1
c
2
 
c     $Id: bead_list.F 19708 2010-10-29 18:04:21Z d3y133 $
 
2
c     $Id: bead_list.F 22987 2012-10-15 16:38:16Z bylaska $
3
3
c
4
4
c
5
5
c     This file contains a list of beads, where a bead is defined
63
63
      return
64
64
      end
65
65
 
 
66
      subroutine set_rtdb_bead_list(rtdb)
 
67
      implicit none
 
68
      integer rtdb
 
69
 
 
70
      integer bead_rtdb
 
71
      common /bead_list/ bead_rtdb
 
72
 
 
73
      bead_rtdb = rtdb
 
74
      return
 
75
      end
 
76
      
 
77
 
 
78
      subroutine reset_bead_list(tag)
 
79
      implicit none
 
80
      integer       rtdb
 
81
      character*(*) tag
 
82
 
 
83
#include "rtdb.fh"
 
84
#include "mafdecls.fh"
 
85
 
 
86
*     ***** local variables ****
 
87
      logical value
 
88
      integer size,taglen
 
89
      character*255 rtdb_name
 
90
 
 
91
      integer bead_rtdb
 
92
      common /bead_list/ bead_rtdb
 
93
 
 
94
*     **** external functions ****
 
95
      integer  inp_strlen
 
96
      external inp_strlen
 
97
 
 
98
      taglen    = inp_strlen(tag)
 
99
      size = 0
 
100
      rtdb_name = tag(1:taglen)//':size'
 
101
      value = rtdb_put(bead_rtdb,rtdb_name,mt_int,1,size) 
 
102
 
 
103
      return
 
104
      end
 
105
 
66
106
c     ***************************************************
67
107
c     *                                                 *
68
108
c     *                add_bead_list                    *
87
127
      logical value
88
128
      integer size,taglen,movecslen,geomlen
89
129
      character*255 rtdb_name
 
130
      real*8 stress(9)
90
131
 
91
132
      integer bead_rtdb
92
133
      common /bead_list/ bead_rtdb
122
163
      value = value.and.
123
164
     >        rtdb_put(bead_rtdb,rtdb_name,mt_dbl,1,0.0d0)
124
165
 
 
166
      call dcopy(9,0.0d0,0,stress,1)
 
167
      rtdb_name = tag(1:taglen)//bead_index_name(size)//':stress'
 
168
      value = value.and.
 
169
     >        rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress)
 
170
 
125
171
      if (.not.value) call errquit('add_bead_list failed',0,0)
126
172
 
127
173
      return
153
199
*     ***** local variables ****
154
200
      logical value
155
201
      integer size,taglen,tmplen,j
156
 
      real*8  energy
 
202
      real*8  energy,stress(9)
157
203
      character*255 rtdb_name,tmp_name
158
204
 
159
205
      integer bead_rtdb
172
218
      value = rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size) 
173
219
      if (i.gt.size) return
174
220
 
 
221
      value = value.and.rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size-1) 
 
222
 
175
223
      do j=i,size-1
176
224
 
177
225
        rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':movecs_name'
178
 
        value = value.and.
179
 
     >          rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)
 
226
        if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name))
 
227
     >     tmp_name = 'bead'//bead_index_name(j+1)//'.movecs'
180
228
        tmplen = inp_strlen(tmp_name)
181
229
        rtdb_name = tag(1:taglen)//bead_index_name(j)//':movecs_name'
182
230
        value = value.and.
183
231
     >          rtdb_cput(bead_rtdb,rtdb_name,1,tmp_name(1:tmplen))
184
232
 
185
233
        rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':geom_name'
186
 
        value = value.and.
187
 
     >          rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name)
 
234
        if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,tmp_name))
 
235
     >     tmp_name = 'bead'//bead_index_name(j+1)//':geom'
188
236
        tmplen = inp_strlen(tmp_name)
189
237
        rtdb_name = tag(1:taglen)//bead_index_name(j)//':geom_name'
190
238
        value = value.and.
197
245
        value = value.and.
198
246
     >          rtdb_put(bead_rtdb,rtdb_name,mt_dbl,1,energy)
199
247
 
 
248
        rtdb_name = tag(1:taglen)//bead_index_name(j+1)//':stress'
 
249
        value = value.and.
 
250
     >          rtdb_get(bead_rtdb,rtdb_name,mt_dbl,9,stress)
 
251
        rtdb_name = tag(1:taglen)//bead_index_name(j)//':stress'
 
252
        value = value.and.
 
253
     >          rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress)
 
254
 
200
255
      end do
201
256
 
202
257
      if (.not.value) call errquit('delete_bead_list failed',0,0)
270
325
#include "geom.fh"
271
326
 
272
327
*     ***** local variables ****
273
 
      logical value
 
328
      logical value,ostress
274
329
      integer size,taglen,nion,geom,geomlen
275
330
      character*255 rtdb_name,geom_name
276
331
 
278
333
      common /bead_list/ bead_rtdb
279
334
 
280
335
*     **** external functions ****
 
336
      logical     bead_includestress
281
337
      integer     inp_strlen,size_bead_list
282
338
      character*7 bead_index_name
 
339
      external    bead_includestress
283
340
      external    inp_strlen,size_bead_list
284
341
      external    bead_index_name
285
342
 
292
349
      end if
293
350
 
294
351
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
295
 
      value = rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
352
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
353
     >   geom_name   = 'bead'//bead_index_name(i)//':geom'
296
354
      geomlen = inp_strlen(geom_name)            
297
355
 
298
 
      value = value.and.geom_create(geom,'bead_tmp') 
 
356
      value =           geom_create(geom,'bead_tmp') 
299
357
      value = value.and.geom_rtdb_load(bead_rtdb,geom,
300
358
     >                                 geom_name(1:geomlen))
301
359
      value = value.and.geom_ncent(geom,nion)
302
360
      value = value.and.geom_destroy(geom)
303
361
      if (.not.value) call errquit('nion_bead_list failed',0,0)
 
362
 
 
363
      if (bead_includestress()) nion = nion + 3
304
364
    
305
365
      nion_bead_list = nion
306
366
      return
328
388
#include "mafdecls.fh"
329
389
#include "geom.fh"
330
390
#include "global.fh"
 
391
#include "errquit.fh"
331
392
 
332
393
*     ***** local variables ****
333
394
      logical value,status
334
395
      integer geom,gradient(2),nion,nelem
335
396
      integer size,taglen,permlen,movecslen,geomlen
336
 
      real*8  energy
 
397
      real*8  energy,stress(9)
337
398
      character*255 rtdb_name,perm_name,movecs_name,geom_name
 
399
      character*32 theory
338
400
 
339
401
      integer bead_rtdb
340
402
      common /bead_list/ bead_rtdb
341
403
 
342
404
*     **** external functions ****
343
 
      logical     task_gradient_gen
344
 
      external    task_gradient_gen
 
405
      logical     task_gradient_gen,bead_includestress
 
406
      external    task_gradient_gen,bead_includestress
345
407
      character*7 bead_index_name
346
408
      integer     inp_strlen
347
409
      external    bead_index_name
351
413
      taglen    = inp_strlen(tag)
352
414
 
353
415
      rtdb_name = tag(1:taglen)//':perm_movecs'
354
 
      value = value.and.
355
 
     >        rtdb_cget(bead_rtdb,rtdb_name,1,perm_name)
 
416
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,perm_name))
 
417
     >   call util_file_prefix('movecs',perm_name)
356
418
      call util_file_name_resolve(perm_name, .false.)
357
419
      permlen = inp_strlen(perm_name)
358
420
 
359
 
      
360
421
      rtdb_name = tag(1:taglen)//':size'
361
 
      value = value.and.
362
 
     >        rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size) 
 
422
      if (.not.rtdb_get(bead_rtdb,rtdb_name,mt_int,1,size))
 
423
     >   size = 0
363
424
      if (i.gt.size) return
364
425
 
365
426
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':movecs_name'
366
 
      value = value.and.
367
 
     >        rtdb_cget(bead_rtdb,rtdb_name,1,movecs_name)
 
427
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,movecs_name))
 
428
     >   movecs_name = 'bead'//bead_index_name(i)//'.movecs'
 
429
 
368
430
      call util_file_name_resolve(movecs_name, .false.)
369
431
      movecslen = inp_strlen(movecs_name)
370
432
 
371
433
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
372
 
      value = value.and.
373
 
     >        rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
434
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
435
     >   geom_name = 'bead'//bead_index_name(i)//':geom'
374
436
      geomlen = inp_strlen(geom_name)
375
437
 
376
438
      value = value.and.
389
451
      end if
390
452
 
391
453
*     *** run gradient task ***
 
454
      if(.not.rtdb_put(bead_rtdb,"scf:converged",mt_log,1,.false.))
 
455
     >  call errquit('scf:converged put',0,0)
 
456
      if(.not.rtdb_put(bead_rtdb,"dft:converged",mt_log,1,.false.))
 
457
     >  call errquit('dft:converged put',0,0)
 
458
 
392
459
      if(.not.rtdb_put(bead_rtdb,"neb:ibead",mt_int,1,i))
393
460
     >  call errquit('neb:ibead put',0,0)
394
 
      if(ga_nodeid().eq.0) 
395
 
     > write(*,*) "neb: running bead ",i
 
461
      if(ga_nodeid().eq.0) write(*,*) "neb: running bead ",i
396
462
      value=value.and.task_gradient_gen(bead_rtdb)
397
463
      if (.not.value) call errquit('run_bead_list failed',1,0)
398
464
 
427
493
      value = value.and.
428
494
     >         rtdb_get(bead_rtdb,'task:gradient',mt_dbl,(3*nion),
429
495
     >                  dbl_mb(gradient(1)))
430
 
     
431
496
 
432
497
      value = value.and.geom_vel_set(geom,dbl_mb(gradient(1)))
433
498
      value = value.and.MA_pop_stack(gradient(2))
435
500
      value = value.and.geom_rtdb_store(bead_rtdb,geom,'geometry')
436
501
      value = value.and.geom_destroy(geom)
437
502
      if (.not.value) call errquit('run_bead_list failed',4,0)
 
503
 
 
504
 
 
505
*     ***** set the stresses ******
 
506
      if (bead_includestress()) then
 
507
 
 
508
         if (.not.rtdb_cget(bead_rtdb, 'task:theory', 1, theory))
 
509
     >   call errquit('run_bead_list: stress theory not specified',
 
510
     >                0,RTDB_ERR)
 
511
         if (theory.eq.'pspw') then
 
512
          if (.not.rtdb_get(bead_rtdb, 'pspw:stress', mt_dbl, 9,stress))
 
513
     >      call errquit('run_bead_list: could not get stress',0,0)
 
514
         else if (theory.eq.'band') then
 
515
          if (.not.rtdb_get(bead_rtdb, 'band:stress', mt_dbl, 9,stress))
 
516
     >       call errquit('run_bead_list: could not get stress',0,0)
 
517
         else if (theory.eq.'paw') then
 
518
          if (.not.rtdb_get(bead_rtdb, 'paw:stress', mt_dbl, 9,stress))
 
519
     >       call errquit('run_bead_list: could not get stress',0,0)
 
520
         else
 
521
           call errquit('run_bead_list: no stress in theory',0,RTDB_ERR)
 
522
         end if
 
523
 
 
524
         rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress'
 
525
         value = value.and.
 
526
     >        rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,stress)
 
527
         if (.not.value) call errquit('run_bead_list failed',2,0)
 
528
      end if
 
529
 
438
530
      if(ga_nodeid().eq.0) then 
439
531
       write(*,*) "neb: finished bead ",i
440
532
       write(*,*) "neb: final energy",energy
594
686
      common /bead_list/ bead_rtdb
595
687
 
596
688
*     **** external functions ****
 
689
      logical     bead_includestress,geom_grad_cart_to_frac
597
690
      character*7 bead_index_name
598
691
      integer     inp_strlen,size_bead_list
 
692
      external    bead_includestress,geom_grad_cart_to_frac
599
693
      external    bead_index_name
600
694
      external    inp_strlen,size_bead_list
601
695
 
605
699
      taglen    = inp_strlen(tag)
606
700
 
607
701
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
608
 
      value     = rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
702
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
703
     >   geom_name = 'bead'//bead_index_name(i)//':geom'
609
704
      geomlen   = inp_strlen(geom_name)            
610
705
 
611
 
      value = value.and.geom_create(geom,'bead_tmp') 
 
706
      value =           geom_create(geom,'bead_tmp') 
612
707
      value = value.and.geom_rtdb_load(bead_rtdb,geom,
613
708
     >                                 geom_name(1:geomlen))
614
709
      value = value.and.geom_ncent(geom,nion)
615
710
      value = value.and.geom_cart_coords_get(geom,c)
 
711
 
 
712
      if (bead_includestress()) then
 
713
 
 
714
*        **** put gradient into fractional ****
 
715
        if (.not. geom_grad_cart_to_frac(geom,c))
 
716
     $  call errquit('coords_get_bead_list: cart_to_frac?',0,0)
 
717
 
 
718
*        **** get stress part of gradient ***
 
719
         if (.not. geom_amatrix_get(geom,c(1,nion+1)))
 
720
     >   call errquit('coords_get_bead_list:failed to get amatrix',0,0)
 
721
      end if
 
722
 
616
723
      value = value.and.geom_destroy(geom)
617
724
      if (.not.value) call errquit('coords_get_bead_list failed',0,0)
618
725
    
654
761
      common /bead_list/ bead_rtdb
655
762
 
656
763
*     **** external functions ****
 
764
      logical     bead_includestress,geom_amatrix_set
657
765
      character*7 bead_index_name
658
766
      integer     inp_strlen,size_bead_list
 
767
      external    bead_includestress,geom_amatrix_set
659
768
      external    bead_index_name
660
769
      external    inp_strlen,size_bead_list
661
770
 
665
774
      taglen    = inp_strlen(tag)
666
775
 
667
776
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
668
 
      value     = rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
777
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
778
     >   geom_name = 'bead'//bead_index_name(i)//':geom'
669
779
      geomlen   = inp_strlen(geom_name)            
670
780
 
671
 
      value = value.and.geom_create(geom,'bead_tmp') 
 
781
      value =           geom_create(geom,'bead_tmp') 
672
782
      value = value.and.geom_rtdb_load(bead_rtdb,geom,
673
783
     >                                 geom_name(1:geomlen))
674
784
      value = value.and.geom_ncent(geom,nion)
 
785
 
 
786
      if (bead_includestress()) then
 
787
 
 
788
        if (.not. geom_frac_to_cart(geom, c))
 
789
     $  call errquit('coords_set_bead_list: frac_to_cart?',0,0)
 
790
 
 
791
        if (.not. geom_amatrix_set(geom,c(1,nion+1)))
 
792
     $  call errquit('coords_set_bead_list:failed to set amatrix',0,0)
 
793
 
 
794
 
 
795
      end if
 
796
 
675
797
      value = value.and.geom_cart_coords_set(geom,c)
676
798
      value = value.and.geom_rtdb_delete(bead_rtdb,geom_name(1:geomlen))
677
799
      value = value.and.geom_rtdb_store(bead_rtdb,geom,
678
800
     >                                  geom_name(1:geomlen))
679
801
      value = value.and.geom_destroy(geom)
680
 
      if (.not.value) call errquit('coords_get_bead_list failed',0,0)
 
802
      if (.not.value) call errquit('coords_set_bead_list failed',0,0)
681
803
    
682
804
      return
683
805
      end
717
839
      common /bead_list/ bead_rtdb
718
840
 
719
841
*     **** external functions ****
 
842
      logical     bead_includestress
720
843
      character*7 bead_index_name
721
844
      integer     inp_strlen,size_bead_list
 
845
      external    bead_includestress
722
846
      external    bead_index_name
723
847
      external    inp_strlen,size_bead_list
724
848
 
728
852
      taglen    = inp_strlen(tag)
729
853
 
730
854
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
731
 
      value     = rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
855
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
856
     >   geom_name   = 'bead'//bead_index_name(i)//':geom'
732
857
      geomlen   = inp_strlen(geom_name)            
733
858
 
734
 
      value = value.and.geom_create(geom,'bead_tmp') 
 
859
      value =           geom_create(geom,'bead_tmp') 
735
860
      value = value.and.geom_rtdb_load(bead_rtdb,geom,
736
861
     >                                 geom_name(1:geomlen))
737
862
      value = value.and.geom_ncent(geom,nion)
738
863
      value = value.and.geom_vel_get(geom,c)
739
864
      value = value.and.geom_destroy(geom)
740
865
      if (.not.value) call errquit('gradient_get_bead_list failed',0,0)
 
866
 
 
867
      if (bead_includestress()) then
 
868
         rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress'
 
869
         value = value.and.
 
870
     >           rtdb_get(bead_rtdb,rtdb_name,mt_dbl,9,c(1,nion+1))
 
871
      end if
741
872
    
742
873
      return
743
874
      end
777
908
      common /bead_list/ bead_rtdb
778
909
 
779
910
*     **** external functions ****
 
911
      logical     bead_includestress
780
912
      character*7 bead_index_name
781
913
      integer     inp_strlen,size_bead_list
 
914
      external    bead_includestress
782
915
      external    bead_index_name
783
916
      external    inp_strlen,size_bead_list
784
917
 
788
921
      taglen    = inp_strlen(tag)
789
922
 
790
923
      rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
791
 
      value     = rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
924
      if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
925
     >   geom_name = 'bead'//bead_index_name(i)//':geom'
792
926
      geomlen   = inp_strlen(geom_name)            
793
927
 
794
 
      value = value.and.geom_create(geom,'bead_tmp') 
 
928
      value =           geom_create(geom,'bead_tmp') 
795
929
      value = value.and.geom_rtdb_load(bead_rtdb,
796
930
     >                                 geom,geom_name(1:geomlen))
797
931
      value = value.and.geom_ncent(geom,nion)
801
935
      value = value.and.geom_rtdb_store(bead_rtdb,
802
936
     >                                  geom,geom_name(1:geomlen))
803
937
      value = value.and.geom_destroy(geom)
 
938
 
 
939
      if (bead_includestress()) then
 
940
         rtdb_name = tag(1:taglen)//bead_index_name(i)//':stress'
 
941
         value = value.and.
 
942
     >           rtdb_put(bead_rtdb,rtdb_name,mt_dbl,9,c(1,nion+1))
 
943
 
 
944
      end if
804
945
      if (.not.value) call errquit('gradient_get_bead_list failed',0,0)
805
946
    
806
947
      return
822
963
c     Entry - tag: name of bead list
823
964
c
824
965
c 
825
 
      subroutine create_xyz_file_bead_list(tag)
 
966
      subroutine create_xyz_file_bead_list(luout,tag,header)
826
967
      implicit none
 
968
      integer luout
827
969
      character*(*) tag
 
970
      logical header
828
971
  
829
972
#include "rtdb.fh"
830
973
#include "mafdecls.fh"
849
992
      external    bead_index_name
850
993
 
851
994
    
852
 
      if (ga_nodeid().eq.0) then
853
 
         write(*,*)
854
 
         write(*,*) 'XYZ FILE for bead_list:',tag
855
 
         write(*,*) '------------------------------------------'
 
995
      if ((ga_nodeid().eq.0).and.(header)) then
 
996
         write(luout,*)
 
997
         write(luout,*) 'XYZ FILE for bead_list:',tag
 
998
         write(luout,*) '------------------------------------------'
856
999
      end if
857
1000
      taglen    = inp_strlen(tag)
858
1001
      nbeads    = size_bead_list(tag)
863
1006
     >                              mt_dbl,1,energy)
864
1007
 
865
1008
         rtdb_name = tag(1:taglen)//bead_index_name(i)//':geom_name'
866
 
         value = value.and.
867
 
     >        rtdb_cget(bead_rtdb,rtdb_name,1,geom_name)
 
1009
         if (.not.rtdb_cget(bead_rtdb,rtdb_name,1,geom_name))
 
1010
     >      geom_name = 'bead'//bead_index_name(i)//':geom'
868
1011
         geomlen = inp_strlen(geom_name)      
869
1012
         value = value.and.geom_rtdb_load(bead_rtdb,
870
1013
     >                             geom,geom_name(1:geomlen))
871
1014
         value = value.and.geom_ncent(geom,nion)
872
1015
 
873
1016
         if (ga_nodeid().eq.0) then
874
 
         write(*,*) nion
875
 
         write(*,*) 'energy=',energy
 
1017
         write(luout,*) nion
 
1018
         write(luout,*) 'energy=',energy
876
1019
         end if
877
1020
         do ii=1,nion       
878
1021
            value = value.and.geom_cent_get(geom,ii,t,rxyz,q)
881
1024
            rxyz(2)= rxyz(2)*0.529177d0
882
1025
            rxyz(3)= rxyz(3)*0.529177d0
883
1026
            if (ga_nodeid().eq.0) then
884
 
            write(*,'(A2,6x,3F12.6)') symbol,rxyz
 
1027
            write(luout,'(A2,6x,3F12.6)') symbol,rxyz
885
1028
            end if
886
1029
         end do
887
1030
      end do
888
1031
      value = value.and.geom_destroy(geom)
889
1032
 
890
 
      if (ga_nodeid().eq.0) then
891
 
         write(*,*) 
 
1033
      if ((ga_nodeid().eq.0).and.(header)) then
 
1034
         write(luout,*) 
892
1035
      end if
893
1036
 
894
1037
      if (.not.value) 
898
1041
      end
899
1042
 
900
1043
 
901
 
 
902
 
 
 
1044
c     ***********************************************
 
1045
c     *                                             *
 
1046
c     *            bead_index_name                  *
 
1047
c     *                                             *
 
1048
c     ***********************************************
903
1049
      character*7 function bead_index_name(i)
904
1050
      integer i
905
1051
 
930
1076
      bead_index_name = name
931
1077
      return
932
1078
      end
 
1079
 
 
1080
c     ***************************************************
 
1081
c     *                                                 *
 
1082
c     *            bead_includestress                   *
 
1083
c     *                                                 *
 
1084
c     ***************************************************
 
1085
      logical function bead_includestress()
 
1086
      implicit none
 
1087
 
 
1088
#include "rtdb.fh"
 
1089
#include "mafdecls.fh"
 
1090
#include "geom.fh"
 
1091
 
 
1092
*     ***** local variables ****
 
1093
      logical ostress
 
1094
 
 
1095
      integer bead_rtdb
 
1096
      common /bead_list/ bead_rtdb
 
1097
 
 
1098
      if (.not.rtdb_get(bead_rtdb,'includestress',mt_log,1,ostress))
 
1099
     >   ostress = .false.
 
1100
 
 
1101
      bead_includestress = ostress
 
1102
      return
 
1103
      end