← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/optim/neb/neb_drv.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/neb_drv.F
1
1
*
2
 
* $Id: neb_drv.F 19708 2010-10-29 18:04:21Z d3y133 $
 
2
* $Id: neb_drv.F 24218 2013-05-15 06:34:07Z bylaska $
3
3
*
4
4
 
 
5
*     ****************************************
 
6
*     *                                      *
 
7
*     *                neb                   *
 
8
*     *                                      *
 
9
*     ****************************************
5
10
      logical function neb(rtdb)
6
11
      implicit none
7
12
      integer rtdb
10
15
#include "mafdecls.fh"
11
16
#include "rtdb.fh"
12
17
#include "util.fh"
 
18
#include "stdio.fh"
 
19
#include "errquit.fh"
13
20
 
14
21
c     
15
22
*     **** local variables ****     
16
 
      logical value, newchain, svalue,verlet,oprint
 
23
      logical value,newchain,svalue,verlet,oprint,done,converged,found
 
24
      logical ostress,finishedstep,converged1
 
25
      character*9 cvg1,cvg2,cvg3,cvg4
17
26
      character*4    mark
18
 
      character*50 bead_list
19
 
      integer ii
20
 
      integer i,it,nbeads,nion,ng,iterNEB,nebsteps
 
27
      character*50 bead_list,filename
 
28
      character*80 title,neb_movecs
 
29
      character*30 tmp
 
30
      character*255 full_filename,full_filename2,filename2
 
31
      character ch_tmp
 
32
      integer ii,shift,print_shift,print_count,jj,manyjj,itm,goodjj
 
33
      integer i,it,nbeads,nion,ng,m,nebsteps,neb_algorithm
 
34
      integer nbeads1,nbeads2
21
35
      integer e1(2),g0(2),g1(2),s(2),t1(2),v1(2),c0(2),c1(2)
22
 
      integer mass(2),dti(2)
23
 
      real*8  path_energy,path_distance,norm,norm0,time_step
24
 
      real*8  Gmax,Grms,Xmax,Xrms,dE,path_energy0
 
36
      integer mass(2),dti(2),cs(2),gs(2)
 
37
      real*8  path_energy,path_distance,norm,norm0,time_step,kbeads
 
38
      real*8  Gmax,Grms,Xmax,Xrms,dE,path_energy0,emid,emax,emin,trust
 
39
      real*8  nebGmax,nebGrms,nebXmax,nebXrms,nebdE,sum,sum2,t,alpha
 
40
      real*8  alphamin,deltaE,sum0,sum0_old,time_step0
 
41
      real*8  emid0,emin0,emax0
25
42
 
26
43
*     **** external functions ****
27
44
      logical task_gradient
30
47
      real*8   energy_bead_list
31
48
      external size_bead_list,nion_bead_list
32
49
      external energy_bead_list
 
50
      character*7 bead_index_name
 
51
      external    bead_index_name
 
52
 
33
53
 
34
54
      oprint = ga_nodeid() .eq. 0
35
55
 
 
56
      if (oprint) then
 
57
         write(luout,*)
 
58
         write(luout,*)
 
59
         call util_print_centered(luout,
 
60
     >        'NWChem Minimum Energy Pathway Program (NEB)',
 
61
     >        40,.true.)
 
62
         write(luout,*)
 
63
         write(luout,*)
 
64
      endif
 
65
 
 
66
      if (rtdb_cget(rtdb,'title',1,title)) then
 
67
         if (oprint) then
 
68
            write(luout,*)
 
69
            write(luout,*)
 
70
            call util_print_centered(6, title, 40, .false.)
 
71
            write(luout,*)
 
72
            write(luout,*)
 
73
         endif
 
74
      endif
 
75
      if (.not. rtdb_get(rtdb,'neb:gmax',mt_dbl,1,nebgmax))
 
76
     >   nebgmax = 0.00045d0
 
77
      if (.not. rtdb_get(rtdb,'neb:grms',mt_dbl,1,nebgrms))
 
78
     >   nebgrms = 0.00030d0
 
79
      if (.not. rtdb_get(rtdb,'neb:xmax',mt_dbl,1,nebxmax))
 
80
     >   nebxmax = 0.000180d0
 
81
      if (.not. rtdb_get(rtdb,'neb:xrms',mt_dbl,1,nebxrms))
 
82
     >   nebxrms = 0.000120d0
 
83
 
36
84
      if (.not.rtdb_get(rtdb,'neb:stepsize',mt_dbl,1,time_step))
37
 
     >  time_step = 10.0d0
 
85
     >  time_step = 1.0d0
 
86
      if (.not.rtdb_get(rtdb,'neb:trust',mt_dbl,1,trust))
 
87
     >  trust = 0.1d0
 
88
      if (.not.rtdb_get(rtdb,'neb:kbeads',mt_dbl,1,kbeads))
 
89
     >  kbeads = 0.1d0
 
90
      if (.not.rtdb_get(rtdb,'neb:steps',mt_int,1,nebsteps))      
 
91
     >   nebsteps = 5
 
92
 
 
93
      if (.not.rtdb_cget(rtdb,'neb:movecs',1,neb_movecs)) then
 
94
         call util_file_prefix('movecs',neb_movecs)
 
95
      end if
 
96
 
 
97
*     *** neb_algorithm ***
 
98
      if (.not.rtdb_get(rtdb,'neb:algorithm',mt_int,1,neb_algorithm))
 
99
     >  neb_algorithm = 0
 
100
      if (.not.rtdb_get(rtdb,'neb:m',mt_int,1,m))  m = 5
 
101
 
 
102
      if (.not.rtdb_get(rtdb,'neb:print_shift',mt_int,1,print_shift))
 
103
     >   print_shift = 0
 
104
 
 
105
*     **** neb_epath filename ****
 
106
      call util_file_prefix('neb_epath',filename)
 
107
      call util_file_name_noprefix(filename,.false.,
 
108
     >                             .false.,
 
109
     >                             full_filename)
 
110
 
 
111
*     **** includestress? ****
 
112
      if (.not. rtdb_get(rtdb,'includestress',mt_log,1,ostress))
 
113
     >   ostress = .false.
 
114
 
38
115
 
39
116
 
40
117
*  RRR only initialize if this is a new neb chain!
41
118
*     **** initialize neb list ****
42
 
      bead_list = 'neb_list'
 
119
      bead_list = 'bead_list'
43
120
      newchain = .false.
44
 
      if (.not.rtdb_get(rtdb,'neb:nebnew',mt_log,1,newchain))
 
121
      if (.not.rtdb_get(rtdb,'bead_list:new',mt_log,1,newchain))
45
122
     >  newchain = .true.
46
123
 
47
124
      if(newchain) then
48
 
        if (ga_nodeid().eq.0)
49
 
     >       write(*,*)'NEW NEB CHAIN, INITIALIZING'
50
 
 
 
125
        if (oprint) write(*,*)'NEW NEB CHAIN, INITIALIZING'
51
126
        call neb_initialize(rtdb,bead_list)
 
127
        print_count = 0
52
128
      else
53
 
         if (ga_nodeid().eq.0)
54
 
     >       write(*,*)'EXISTING NEB CHAIN? RESTARTING'
55
 
 
56
 
      endif
 
129
         if (oprint) write(*,*)'EXISTING NEB CHAIN? RESTARTING'
 
130
         call set_rtdb_bead_list(rtdb)
 
131
         nbeads1 = size_bead_list(bead_list)
 
132
         if (.not.rtdb_get(rtdb,'neb:nbeads',mt_int,1,nbeads2))
 
133
     >      nbeads2 = 5
 
134
         if (nbeads1.ne.nbeads2) then
 
135
            if (oprint) write(*,*) 'RESIZING NEB CHAIN'
 
136
            call neb_resize_path(rtdb,bead_list,nbeads1,nbeads2)
 
137
         end if
 
138
         if (.not.rtdb_get(rtdb,'neb:print_count',mt_int,1,print_count))
 
139
     >      print_count = 0
 
140
      end if
57
141
      newchain = .false.
58
142
 
59
143
      nbeads = size_bead_list(bead_list)
60
144
      nion   = nion_bead_list(bead_list,1)
61
145
      ng     = 3*nion*nbeads
62
146
 
63
 
*     *** is verlet algorithm used ***
64
 
      if (.not.rtdb_get(rtdb,'neb:verlet',mt_log,1,verlet))
65
 
     >  verlet = .false.
 
147
      if (oprint) then
 
148
 
 
149
         if (neb_algorithm.eq.0) then
 
150
            tmp = "Fixed Point"
 
151
         else if (neb_algorithm.eq.1) then
 
152
            tmp = "Verlet"
 
153
         else if (neb_algorithm.eq.2) then
 
154
            tmp = "refining conjugate gradient"
 
155
         else
 
156
            tmp = "not implemented"
 
157
         end if
 
158
         
 
159
         write(luout,*)
 
160
         write(luout,1) nebgmax,nebgrms,nebxmax,nebxrms,
 
161
     >                  time_step,trust,nebsteps,nbeads,m,nion,kbeads,
 
162
     >                  neb_algorithm,tmp,
 
163
     >                  neb_movecs
 
164
         if (ostress) 
 
165
     >      write(luout,*) ' INCLUDING STRESS !!!!!!!!!!!!!!!!'
 
166
         write(luout,*)
 
167
         write(luout,*)
 
168
 1       format(
 
169
     >        ' maximum gradient threshold         (gmax) = ', f10.6,/,
 
170
     >        ' rms gradient threshold             (grms) = ', f10.6,/,
 
171
     >        ' maximum cartesian step threshold   (xmax) = ', f10.6,/,
 
172
     >        ' rms cartesian step threshold       (xrms) = ', f10.6,/,
 
173
     >        0p,/,
 
174
     >        ' step size                      (stepsize) = ', f10.6,/,
 
175
     >        ' fixed trust radius                (trust) = ', f10.6,/,
 
176
     >        ' maximum number of steps         (maxiter) = ', i4,/,
 
177
     >        ' number of images in path         (nbeads) = ', i4,/,
 
178
     >        ' number of histories               (nhist) = ', i4,/,
 
179
     >        ' number of atoms                           = ', i4,/,
 
180
     >        ' NEB spring constant in a.u.      (kbeads) = ', f10.6,/,
 
181
     >        ' NEB algorithm                 (algorithm) = ', i4,
 
182
     >                                                   1x,'(',a30,')'/,
 
183
     >        ' NEB movecs filename                       = ', a)
 
184
      end if
66
185
 
67
186
 
68
187
*     **** allocate space for gradients and coordinates ****
84
203
      value  = value.and.
85
204
     >         MA_alloc_get(mt_dbl,ng,'dti',dti(2),dti(1))
86
205
      if (.not.value) 
87
 
     >  call errquit('neb failed - increase stack memory',1,0)
 
206
     >  call errquit('neb failed - increase heap memory',1,0)
 
207
 
 
208
 
 
209
      value  = 
 
210
     >         MA_alloc_get(mt_dbl,m*ng,'cs',cs(2),cs(1))
 
211
      value  = value.and.
 
212
     >         MA_alloc_get(mt_dbl,m*ng,'gs',gs(2),gs(1))
 
213
      if (.not.value) 
 
214
     >  call errquit('neb failed - increase heap memory',1,0)
 
215
 
88
216
 
89
217
*     *** set dti ***
90
218
      value  = MA_alloc_get(mt_dbl,nion,'mass',mass(2),mass(1))      
91
219
      if (.not.value) 
92
 
     >  call errquit('neb failed - increase stack memory',2,0)
 
220
     >  call errquit('neb failed - increase heap memory',2,0)
93
221
      call neb_masses_get(rtdb,dbl_mb(mass(1)))
94
222
 
95
223
      do i=1,nbeads
107
235
  
108
236
 
109
237
*     **** initial step ****
110
 
      if (ga_nodeid().eq.0) 
111
 
     >   write(*,*) "neb: Calculating Initial Path Energy "
 
238
      if (oprint) write(luout,*) "neb: Calculating Initial Path Energy"
 
239
      call neb_coords_get(bead_list,dbl_mb(c1(1)))
112
240
      call runall_bead_list(bead_list,task_gradient)
113
241
      call neb_energies_get(bead_list,dbl_mb(e1(1)))
114
 
      call neb_coords_get(bead_list,dbl_mb(c1(1)))
115
 
      call neb_gradient_get(bead_list,
 
242
      call dcopy(ng,0.0d0,0,dbl_mb(t1(1)),1)
 
243
      call dcopy(ng,0.0d0,0,dbl_mb(g1(1)),1)
 
244
      call neb_gradient_get(bead_list,kbeads,
116
245
     >                      dbl_mb(c1(1)),
117
246
     >                      dbl_mb(e1(1)),
118
247
     >                      dbl_mb(t1(1)),
119
248
     >                      dbl_mb(g1(1))) 
 
249
      sum0 = ddot(ng,dbl_mb(g1(1)),1,dbl_mb(g1(1)),1)
 
250
      sum0_old = sum0
 
251
      manyjj = 0
 
252
      goodjj = 0
 
253
      time_step0 = time_step
 
254
 
120
255
 
121
256
      call neb_path_energy(bead_list,
122
257
     >                     path_distance,
123
258
     >                     path_energy) 
124
259
 
125
 
      IF (ga_nodeid().eq.0) THEN
126
 
      write(*,*)
127
 
      write(*,*) "neb: Initial Path Energy    "
128
 
      write(*,*) "neb: -----------------------"
129
 
      do i=1,nbeads
130
 
         write(*,*) "neb: ",i,dbl_mb(e1(1)+i-1)
131
 
      end do
132
 
      write(*,*)
133
 
      END IF
134
 
      call create_xyz_file_bead_list(bead_list)
 
260
      if (oprint) then
 
261
         write(luout,*) "neb: sum0=",sum0,ng
 
262
         write(luout,*)
 
263
         write(luout,*) "neb: Initial Path Energy    "
 
264
         write(luout,*) "neb: -----------------------"
 
265
         do i=1,nbeads
 
266
            write(luout,*) "neb: ",i,dbl_mb(e1(1)+i-1)
 
267
         end do
 
268
         write(luout,*)
 
269
      end if
 
270
      call create_xyz_file_bead_list(luout,bead_list,.true.)
135
271
 
136
272
      norm = dsqrt(ddot(ng,dbl_mb(g1(1)),1,dbl_mb(g1(1)),1))
137
 
      IF (ga_nodeid().eq.0)
138
 
     >       write(*,*) "Path Energy, Path Distance, |G_neb|:",
 
273
      if (oprint) write(luout,*) "Path Energy, Path Distance, |G_neb|:",
139
274
     >            path_energy,path_distance,norm
140
275
 
141
 
 
142
 
 
143
 
      svalue =  rtdb_get(rtdb,'neb:steps',mt_int,1,nebsteps)      
144
 
      if(svalue)then
145
 
      iterNEB = nebsteps
146
 
      IF (ga_nodeid().eq.0)
147
 
     >      write(*,*)'NEB iterations  =',iterNEB
148
 
      else
149
 
         iterNEB = 5
150
 
      IF (ga_nodeid().eq.0)
151
 
     >      write(*,*)'SHORTRUN, NEB iterations =',iterNEB
152
 
      endif
153
 
 
154
276
      call dcopy(ng,dbl_mb(g1(1)),1,dbl_mb(s(1)),1)
155
277
 
156
 
      do it=1,iterNEB
 
278
      it  = 0
 
279
      itm = 0
 
280
      done = .false.
 
281
      do while (.not.done)
 
282
         it  = it + 1
 
283
         itm = itm + 1
 
284
         if (oprint) write(luout,*) 'neb: iteration #',it,itm
157
285
 
158
 
         IF (ga_nodeid().eq.0)
159
 
     >         write(*,*)'neb: iteration #',it
160
286
*        *** save old forces  and coordinates ***
161
287
         call dcopy(ng,dbl_mb(c1(1)),1,dbl_mb(c0(1)),1)
162
288
         call dcopy(ng,dbl_mb(g1(1)),1,dbl_mb(g0(1)),1)
163
289
 
164
290
         norm0=norm
165
291
 
166
 
          if(verlet) THEN
167
 
            IF (ga_nodeid().eq.0)
168
 
     >         write(*,*)'neb: using verlet algroithm'
 
292
         
 
293
*         ***** fixed point iteration ****
 
294
          if (neb_algorithm.eq.0) then
 
295
             if (oprint) write(luout,*)'neb: using fixed point'
 
296
 
 
297
             if (itm.le.m) then
 
298
                shift = (itm-1)*ng
 
299
                call dcopy(ng,dbl_mb(c0(1)),1,dbl_mb(cs(1)+shift),1)
 
300
                call dcopy(ng,dbl_mb(g0(1)),1,dbl_mb(gs(1)+shift),1)
 
301
                if (itm.gt.1) then
 
302
                   call neb_lmbfgs(ng,itm,
 
303
     >                             dbl_mb(cs(1)),
 
304
     >                             dbl_mb(gs(1)),
 
305
     >                             dbl_mb(s(1)))
 
306
                else
 
307
                   call dcopy(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
 
308
                end if
 
309
             else
 
310
                do ii=1,m-1
 
311
                   shift = (ii-1)*ng
 
312
                   call dcopy(ng,dbl_mb(cs(1)+shift+ng),1,
 
313
     >                           dbl_mb(cs(1)+shift),   1)
 
314
                   call dcopy(ng,dbl_mb(gs(1)+shift+ng),1,
 
315
     >                           dbl_mb(gs(1)+shift),   1)
 
316
 
 
317
                end do
 
318
                shift = (m-1)*ng
 
319
                call dcopy(ng,dbl_mb(c0(1)),1,dbl_mb(cs(1)+shift),1)
 
320
                call dcopy(ng,dbl_mb(g0(1)),1,dbl_mb(gs(1)+shift),1)
 
321
                call neb_lmbfgs(ng,m,
 
322
     >                          dbl_mb(cs(1)),
 
323
     >                          dbl_mb(gs(1)),
 
324
     >                          dbl_mb(s(1)))
 
325
             end if
 
326
c             call neb_project_gradient(nion,nbeads,
 
327
c     >                                dbl_mb(t1(1)),
 
328
c     >                                dbl_mb(s(1)))
 
329
 
 
330
             finishedstep = .false.
 
331
             alpha = time_step
 
332
             jj = 0
 
333
             do while ((.not.finishedstep).and.(jj.lt.9))
 
334
             jj = jj + 1
 
335
             sum = ddot(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
 
336
 
 
337
             if ((sum.le.0.0d0).or.(sum0.gt.sum0_old))  then
 
338
                call dcopy(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
 
339
                itm = 0
 
340
                if (oprint) write(*,*) "neb: s=g and itm reset to 0"
 
341
             end if
 
342
 
 
343
             sum2 = dsqrt(ddot(ng,dbl_mb(s(1)),1,dbl_mb(s(1)),1))
 
344
     >              /dble(nbeads)
 
345
             if (oprint) 
 
346
     >          write(*,*) "neb: |<s|s>|,<g|s>=",sum2,sum
 
347
 
 
348
             call neb_move(ng,
 
349
     >                     (-alpha),
 
350
     >                     dbl_mb(c0(1)),
 
351
     >                     dbl_mb(c1(1)),
 
352
     >                     dbl_mb(s(1)))
 
353
              call neb_coords_set(bead_list,dbl_mb(c1(1)))
 
354
 
 
355
              if (oprint) write(luout,*) "neb: running internal beads"
 
356
              call runmid_bead_list(bead_list,task_gradient)
 
357
 
 
358
              call neb_energies_get(bead_list,dbl_mb(e1(1)))
 
359
              call neb_gradient_get(bead_list,kbeads,
 
360
     >                             dbl_mb(c1(1)),
 
361
     >                             dbl_mb(e1(1)),
 
362
     >                             dbl_mb(t1(1)),
 
363
     >                             dbl_mb(g1(1)))
 
364
               sum0     = ddot(ng,dbl_mb(g1(1)),1,dbl_mb(g1(1)),1)
 
365
 
 
366
               finishedstep = (sum0.le.sum0_old)
 
367
               if (.not.finishedstep) then
 
368
                   alpha = 0.50d0*alpha
 
369
                   manyjj = manyjj+1
 
370
                   goodjj = 0
 
371
                   if (manyjj.gt.2) then
 
372
                      time_step = 0.5d0*time_step
 
373
                      manyjj = 0
 
374
                      if (oprint) 
 
375
     >                   write(*,*) "neb: reducing timestep=",time_step
 
376
                   end if
 
377
               end if
 
378
               if (finishedstep.and.(jj.eq.1)) then
 
379
                  goodjj = goodjj + 1
 
380
                  if (goodjj.gt.10) then
 
381
                     time_step = 2.0d0*time_step
 
382
                     if (time_step.gt.time_step0) time_step = time_step0
 
383
                     goodjj = 0
 
384
                  end if
 
385
               end if
 
386
 
 
387
               if(oprint) 
 
388
     >             write(*,*) "neb: sum0,sum0_old=",
 
389
     >                        sum0,sum0_old,jj,finishedstep,alpha
 
390
 
 
391
c               if (finishedstep) sum0_old = sum0
 
392
 
 
393
             end do
 
394
             sum0_old = sum0
 
395
 
 
396
 
 
397
*         ***** Verlet iteration ****
 
398
          else if (neb_algorithm.eq.1) then
 
399
            IF (oprint) write(luout,*)'neb: using verlet algroithm'
169
400
            call neb_verlet_update(ng,
170
401
     >                         dbl_mb(c0(1)),
171
402
     >                         dbl_mb(c1(1)),
176
407
            call neb_coords_set(bead_list,dbl_mb(c1(1)))
177
408
            call runmid_bead_list(bead_list,task_gradient)
178
409
            call neb_energies_get(bead_list,dbl_mb(e1(1)))
179
 
            call neb_gradient_get(bead_list,
 
410
            call neb_gradient_get(bead_list,kbeads,
180
411
     >                            dbl_mb(c1(1)),
181
412
     >                            dbl_mb(e1(1)),
182
413
     >                            dbl_mb(t1(1)),
183
414
     >                            dbl_mb(g1(1)))
184
415
 
185
 
          else 
186
 
            if(oprint) write(*,*)'neb: using cg algorithm'
 
416
*         ***** refined CG iteration ****
 
417
          else if (neb_algorithm.eq.2) then
 
418
            if (oprint) write(luout,*) 'neb: using cg algorithm'
187
419
            call neb_cg_direction(ng,
188
420
     >                         dbl_mb(g0(1)),
189
421
     >                         dbl_mb(g1(1)),
191
423
 
192
424
            do ii=1,10
193
425
 
194
 
            if(oprint)  write(*,*) "neb: refining time step"
 
426
            if (oprint) write(luout,*) "neb: refining time step"
195
427
              call neb_move(ng,
196
428
     >                      time_step,
197
429
     >                      dbl_mb(c0(1)),
199
431
     >                      dbl_mb(s(1)))
200
432
 
201
433
              call neb_coords_set(bead_list,dbl_mb(c1(1)))
202
 
              if(oprint)  write(*,*) "neb: running internal beads"
 
434
              if (oprint) write(luout,*) "neb: running internal beads"
203
435
              call runmid_bead_list(bead_list,task_gradient)
204
436
              call neb_energies_get(bead_list,dbl_mb(e1(1)))
205
 
              call neb_gradient_get(bead_list,
 
437
              call neb_gradient_get(bead_list,kbeads,
206
438
     >                             dbl_mb(c1(1)),
207
439
     >                             dbl_mb(e1(1)),
208
440
     >                             dbl_mb(t1(1)),
209
441
     >                             dbl_mb(g1(1)))
210
442
 
211
443
              norm = dsqrt(ddot(ng,dbl_mb(g1(1)),1,dbl_mb(g1(1)),1))
212
 
              if(oprint)  write(*,*) "neb: new gnorm=",norm
213
 
              if(oprint)  write(*,*) "neb: old gnorm0=",norm0
 
444
              if(oprint)  write(luout,*) "neb: new gnorm=",norm
 
445
              if(oprint)  write(luout,*) "neb: old gnorm0=",norm0
214
446
              if(norm.gt.norm0) then
215
447
                time_step=time_step/2.0d0
216
 
                if(oprint)  
217
 
     >            write(*,*) "neb: reducing time step",time_step 
 
448
                if (oprint) 
 
449
     >             write(luout,*) "neb: reducing time step ",time_step
218
450
              else
219
451
                call dscal(ng,time_step,dbl_mb(s(1)),1)
220
 
                if(oprint)  
221
 
     >            write(*,*) "neb: accepting time step"
 
452
                if (oprint) 
 
453
     >             write(luout,*) "neb: accepting time step ",time_step
222
454
                go to 19
223
455
              end if
224
 
           
225
456
            end do
226
457
         end if
227
 
 
228
458
19       continue         
229
459
 
 
460
         if (oprint) then
 
461
            inquire(file=full_filename,exist=found)
 
462
*           **** CIF FILE already exists - parse to EOF ****
 
463
            if (found) then
 
464
              open(unit=19,file=full_filename,form='formatted',
 
465
     >             status='old')
 
466
              do while(.true.)
 
467
                read(19,*,ERR=30,END=30) ch_tmp
 
468
              end do
 
469
 30           continue
 
470
#if defined(FUJITSU_SOLARIS) || defined(PSCALE) || defined(SOLARIS) || defined(__crayx1) || defined(GCC46)
 
471
              backspace 19
 
472
#endif
 
473
            write(19,*) " "
 
474
*           **** .neb_epath FILE does not exist ****
 
475
            else
 
476
              open(unit=19,file=full_filename,form='formatted')
 
477
            end if
 
478
            write(19,*) 
 
479
     > "#-------------------------------------------------------"
 
480
            write(19,*)  "# NEB Path iteration = ",it
 
481
            write(19,*)  "# algorithm          = ",neb_algorithm
 
482
            write(19,*)  "# nbeads             = ",nbeads
 
483
            write(19,*)  "# nhist              = ",m
 
484
            write(19,*)  "# natoms             = ",nion
 
485
            write(19,*)  "# kbeads             = ",kbeads
 
486
            write(19,*)  "# stepsize           = ",time_step
 
487
            write(19,*)  "# trust              = ",trust
 
488
            write(19,*)  "# path energy        = ",path_energy
 
489
            write(19,*)  "# path distance      = ",path_distance
 
490
            write(19,*)  "# dE                 = ",dE
 
491
            write(19,*)  "# Gmax               = ",Gmax
 
492
            write(19,*)  "# Grms               = ",Grms
 
493
            write(19,*)  "# Xmax               = ",Xmax
 
494
            write(19,*)  "# Xrms               = ",Xrms
 
495
            write(19,*) 
 
496
     > "#-------------------------------------------------------"
 
497
            do i=1,nbeads
 
498
               t = (i-1)/dble(nbeads-1)
 
499
               write(19,*) t,dbl_mb(e1(1)+i-1)
 
500
            end do
 
501
            close(19)
 
502
         end if
 
503
 
230
504
 
231
505
*        *** RRR write out cumulative path energy
232
 
         IF (ga_nodeid().eq.0) THEN
233
 
         write(*,*)
234
 
         write(*,*) "neb: Path Energy #",it
235
 
         write(*,*) "----------------------------"
236
 
         do i=1,nbeads
237
 
            write(*,*) "neb: ",i,dbl_mb(e1(1)+i-1)
238
 
         end do
239
 
         write(*,*)
240
 
         END IF
241
 
         call create_xyz_file_bead_list(bead_list)
242
 
 
243
 
         path_energy0 = path_energy
 
506
         if (oprint) then
 
507
            write(luout,*)
 
508
            write(luout,*) "neb: Path Energy #",it
 
509
            write(luout,*) "----------------------------"
 
510
            do i=1,nbeads
 
511
               write(luout,*) "neb: ",i,dbl_mb(e1(1)+i-1)
 
512
            end do
 
513
            write(luout,*)
 
514
         end if
 
515
         call create_xyz_file_bead_list(luout,bead_list,.true.)
 
516
 
 
517
*        **** write out current xyz file ***
 
518
         print_count = print_count + 1
 
519
         if (print_shift.gt.0) then
 
520
            if (mod(print_count,print_shift).eq.0) then
 
521
               call util_file_prefix(
 
522
     >                  'nebpath'//bead_index_name(print_count)//'.xyz',
 
523
     >                  filename2)
 
524
               call util_file_name_noprefix(filename2,.false.,
 
525
     >                             .false.,
 
526
     >                             full_filename2)
 
527
                if (oprint) 
 
528
     >            open(unit=23,file=full_filename2,form='formatted')
 
529
                call create_xyz_file_bead_list(23,bead_list,.false.)
 
530
                if (oprint) close(23)
 
531
            end if
 
532
         end if
 
533
         if (.not.rtdb_put(rtdb,'neb:print_count',mt_int,1,print_count))
 
534
     >      call errquit('setting neb:print_count failed',4,RTDB_ERR)
 
535
 
 
536
         emid0 = path_energy/path_distance
244
537
         call neb_path_energy(bead_list,
245
538
     >                        path_distance,
246
539
     >                        path_energy) 
247
 
         dE = path_energy - path_energy0
 
540
         dE = dabs((path_energy/path_distance) - emid0)
248
541
         call neb_calc_convergence(ng,dbl_mb(g1(1)),
249
542
     >                                dbl_mb(c0(1)),
250
543
     >                                dbl_mb(c1(1)),
251
544
     >                                Gmax,Grms,Xmax,Xrms)
252
545
 
253
 
!        IF (ga_nodeid().eq.0) THEN
254
 
!             write(*,*) "@  Iteration#:",it
255
 
!             write(*,*) "@  Path Energy:",path_energy
256
 
!             write(*,*) "@  Path Distance:",path_distance
257
 
!             write(*,*) "@  |G_neb|:",norm
258
 
!        END IF
 
546
         emid0 = emid
 
547
         emin0 = emin
 
548
         emax0 = emax
 
549
 
 
550
         emin = +99.0d99
 
551
         emax = -99.0d99
 
552
         do i=1,nbeads
 
553
             if (dbl_mb(e1(1)+i-1)>emax) 
 
554
     >          emax = dbl_mb(e1(1)+i-1)
 
555
             if (dbl_mb(e1(1)+i-1)<emin) 
 
556
     >          emin = dbl_mb(e1(1)+i-1)
 
557
         end do
 
558
         i = nbeads/2
 
559
         if (i.lt.1) i = 1
 
560
         emid = dbl_mb(e1(1)+i-1)
 
561
 
 
562
         dE = dE + dabs(emid-emid0)
 
563
         dE = dE + dabs(emin-emin0)
 
564
         dE = dE + dabs(emax-emax0)
259
565
 
260
566
         if (oprint) then
 
567
           cvg1 = ' '
 
568
           cvg2 = ' '
 
569
           cvg3 = ' '
 
570
           cvg4 = ' '
 
571
           if (gmax .lt. nebGmax) cvg1 = '     ok  '
 
572
           if (grms .lt. nebGrms) cvg2 = '     ok  '
 
573
           if (xrms .lt. nebXrms) cvg3 = '     ok  '
 
574
           if (xmax .lt. nebXmax) cvg4 = '     ok  '
 
575
 
261
576
           mark = '@neb'
262
577
           if (it .gt. 1) mark = ' '
263
 
           write(6,1) mark, mark
 
578
 
 
579
           write(luout,12)  mark," "
 
580
           write(luout,12)  mark,"NEB Method"
 
581
           write(luout,13)  mark,"algorithm      = ",neb_algorithm
 
582
           write(luout,13)  mark,"maxiter        = ",nebsteps
 
583
           write(luout,13)  mark,"nbeads         = ",nbeads
 
584
           write(luout,13)  mark,"nhist          = ",m
 
585
           write(luout,13)  mark,"natoms         = ",nion
 
586
           write(luout,14)  mark,"stepsize       = ",time_step
 
587
           write(luout,14)  mark,"trust          = ",trust
 
588
           write(luout,14)  mark,"kbeads         = ",kbeads
 
589
           write(luout,14)  mark,"Gmax tolerance = ",nebGmax
 
590
           write(luout,14)  mark,"Grms tolerance = ",nebGrms
 
591
           write(luout,14)  mark,"Xmax tolerance = ",nebXmax
 
592
           write(luout,14)  mark,"Xrms tolerance = ",nebXrms
 
593
           write(luout,12)  mark," "
 
594
 
 
595
           write(luout,10) mark, mark
 
596
     
264
597
           mark = '@neb'
265
 
           write(6,2) mark, it-1, path_energy, dE,
266
 
     $       Gmax, Grms, Xrms, Xmax, util_wallsec()
267
 
 1         format(
268
 
     $        /,a4,' Step     Path Energy   Delta E   Gmax',
269
 
     $        '     Grms     Xrms     Xmax   Walltime',
270
 
     $        /,a4,' ---- ---------------- -------- --------',
271
 
     $        ' -------- -------- -------- --------')
272
 
 2         format(
273
 
     $        a4,i5,f17.8,1p,d9.1,0p,4f9.5,f9.1,/,
274
 
     $        1x,5x,17x,9x,4a9,/)
 
598
           write(luout,11) mark, it, path_energy/path_distance, 
 
599
     >       emid,emin,emax,
 
600
     >       Gmax, Grms, Xrms, Xmax, util_wallsec(),
 
601
     >       cvg1,cvg2,cvg3,cvg4
 
602
 10        format(
 
603
     >        /,a4,' Step    Intrinsic E    Mid-Point E ',
 
604
     >        '     Minimum E      Maximum E   Gmax',
 
605
     >        '     Grms     Xrms     Xmax   Walltime',
 
606
     >        /,a4,' ---- -------------- -------------- --------------',
 
607
     >        ' --------------',
 
608
     >        ' -------- -------- -------- -------- --------')
 
609
 11        format(
 
610
     >        a4,i5,4f15.6,1p,4f9.5,f9.1,/,
 
611
     >        1x,5x,17x,9x,4a9,/)
 
612
 12        format(a4,1x,a)
 
613
 13        format(a4,1x,a,i9)
 
614
 14        format(a4,1x,a,e9.3)
 
615
 
275
616
         endif
276
617
 
 
618
        converged1 = (Gmax.le.nebgmax)
 
619
        converged =               (Grms.le.nebgrms)
 
620
        converged = converged.and.(Xmax.le.nebxmax)
 
621
        converged = converged.and.(Xrms.le.nebxrms)
 
622
        if (dE.gt.1.0d-6) converged = converged.and.converged1
 
623
 
 
624
        done = converged.or.(it.ge.nebsteps)
277
625
  
278
626
      end do
279
 
 
 
627
      if (oprint) then
 
628
         if (converged) then
 
629
            write(luout,'(a)') "@neb  NEB calculation converged"
 
630
            if (.not.converged1) 
 
631
     >         write(luout,'(2a)') 
 
632
     >           "@neb  However NEB Gmax not converged",
 
633
     >           "...Try increasing the number of beads."
 
634
         else
 
635
            write(luout,'(a)') "@neb   NEB calculation not converged"
 
636
         end if
 
637
          
 
638
      end if
 
639
 
 
640
 
 
641
*     **** write out final path energies ****
 
642
      call util_file_prefix('neb_final_epath',filename)
 
643
      call util_file_name_noprefix(filename,.false.,
 
644
     >                             .false.,
 
645
     >                             full_filename)
 
646
      if (oprint) then
 
647
         open(unit=19,file=full_filename,form='formatted')
 
648
         write(19,*) 
 
649
     > "#-------------------------------------------------------"
 
650
         write(19,*)  "# NEB Path"
 
651
         write(19,*)  "# algorithm     = ",neb_algorithm
 
652
         write(19,*)  "# nbeads        = ",nbeads
 
653
         write(19,*)  "# nhist         = ",m
 
654
         write(19,*)  "# natoms        = ",nion
 
655
         write(19,*)  "# kbeads        = ",kbeads
 
656
         write(19,*)  "# stepsize      = ",time_step
 
657
         write(19,*)  "# trust         = ",trust
 
658
         write(19,*)  "# path energy   = ",path_energy
 
659
         write(19,*)  "# path distance = ",path_distance
 
660
         write(19,*)  "# dE            = ",dE
 
661
         write(19,*)  "# Gmax          = ",Gmax
 
662
         write(19,*)  "# Grms          = ",Grms
 
663
         write(19,*)  "# Xmax          = ",Xmax
 
664
         write(19,*)  "# Xrms          = ",Xrms
 
665
         write(19,*) 
 
666
     > "#-------------------------------------------------------"
 
667
         do i=1,nbeads
 
668
            t = (i-1)/dble(nbeads-1)
 
669
            write(19,*) t,dbl_mb(e1(1)+i-1)
 
670
         end do
 
671
         close(19)
 
672
      end if
 
673
 
 
674
*     **** write out final xyz movies energies ****
 
675
      call util_file_prefix('neb_final.xyz',filename)
 
676
      call util_file_name_noprefix(filename,.false.,
 
677
     >                             .false.,
 
678
     >                             full_filename)
 
679
       if (oprint) open(unit=19,file=full_filename,form='formatted')
 
680
       call create_xyz_file_bead_list(19,bead_list,.false.)
 
681
       if (oprint) close(19)
 
682
 
 
683
 
 
684
      value = value.and.MA_free_heap(cs(2))
 
685
      value = value.and.MA_free_heap(gs(2))
280
686
      value = value.and.MA_free_heap(dti(2))
281
687
      value = value.and.MA_free_heap(c1(2))
282
688
      value = value.and.MA_free_heap(c0(2))
288
694
      value = value.and.MA_free_heap(e1(2))      
289
695
      if (.not.value) call errquit('neb failed',4,0)
290
696
 
 
697
      if (.not.rtdb_put(rtdb,'bead_list:new',mt_log,1,.false.))
 
698
     > call errquit('setting bead_list:new failed',4,RTDB_ERR)
 
699
 
291
700
      call ga_sync()
292
701
      neb = .true. 
293
702
      end