← 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_utils.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_utils.F
1
1
*
2
 
* $Id: neb_utils.F 19708 2010-10-29 18:04:21Z d3y133 $
 
2
* $Id: neb_utils.F 24218 2013-05-15 06:34:07Z bylaska $
3
3
*
 
4
 
 
5
*     ***********************************************
 
6
*     *                                             *
 
7
*     *             neb_path_energy                 *
 
8
*     *                                             *
 
9
*     ***********************************************
4
10
      subroutine neb_path_energy(bead_list,dpath,epath)
5
11
      implicit none
6
12
      character*(*) bead_list
61
67
      return
62
68
      end
63
69
 
64
 
      subroutine neb_gradient_get(bead_list,c,e,t,g)
 
70
*     ***********************************************
 
71
*     *                                             *
 
72
*     *             neb_gradient_get                *
 
73
*     *                                             *
 
74
*     ***********************************************
 
75
      subroutine neb_gradient_get(bead_list,kbeads,c,e,t,g)
65
76
      implicit none
66
77
      character*(*) bead_list
 
78
      real*8 kbeads
67
79
      real*8 c(*)
68
80
      real*8 e(*)
69
81
      real*8 t(*)
80
92
      external ddot
81
93
 
82
94
 
83
 
      k = 0.1d0 !*** what's a good value?
 
95
      k = kbeads
84
96
 
85
97
      nbeads = size_bead_list(bead_list)
86
98
      nion   = nion_bead_list(bead_list,1)
105
117
      return
106
118
      end
107
119
 
 
120
      subroutine neb_gradient_get0(bead_list,g)
 
121
      implicit none
 
122
      character*(*) bead_list
 
123
      real*8 g(*)
 
124
 
 
125
*     **** local variables ****
 
126
      integer i,index,nbeads,nion
 
127
 
 
128
*     **** external functions ****
 
129
      integer  size_bead_list,nion_bead_list
 
130
      external size_bead_list,nion_bead_list
 
131
 
 
132
 
 
133
      nbeads = size_bead_list(bead_list)
 
134
      nion   = nion_bead_list(bead_list,1)
 
135
 
 
136
      call dcopy(nbeads*3*nion,0.0d0,0,g,1)
 
137
      do i=2,(nbeads-1)
 
138
         index = (i-1)*3*nion + 1
 
139
         call gradient_get_bead_list(bead_list,i,g(index))
 
140
      end do
 
141
      return
 
142
      end
 
143
 
 
144
 
 
145
      subroutine neb_project_gradient(nion,nbeads,t,g)
 
146
      implicit none
 
147
      integer nion,nbeads
 
148
      real*8 t(*),g(*)
 
149
 
 
150
*     **** local variables ****
 
151
      integer i,index
 
152
      real*8 norm
 
153
 
 
154
*     **** external functions ****
 
155
      real*8   ddot
 
156
      external ddot
 
157
 
 
158
      do i=2,(nbeads-1)
 
159
         index = (i-1)*3*nion+1
 
160
         norm = ddot(3*nion,t(index),1,g(index),1)
 
161
         call daxpy(3*nion,(-norm),
 
162
     >              t(index),1,
 
163
     >              g(index),1)
 
164
      end do
 
165
 
 
166
      return
 
167
      end
 
168
     
 
169
  
 
170
 
 
171
*     ***********************************************
 
172
*     *                                             *
 
173
*     *             neb_add_spring_gradient         *
 
174
*     *                                             *
 
175
*     ***********************************************
108
176
      subroutine neb_add_spring_gradient(nbeads,nion,c,t,k,gs)
109
177
      implicit none
110
178
      integer nbeads,nion
151
219
      return
152
220
      end
153
221
 
 
222
*     ***********************************************
 
223
*     *                                             *
 
224
*     *             neb_tangent                     *
 
225
*     *                                             *
 
226
*     ***********************************************
154
227
      subroutine neb_tangent(nbeads,nion,c,e,t)
155
228
      implicit none
156
229
      integer nbeads,nion
183
256
     >               c(index),1,
184
257
     >               t(rp),   1)
185
258
 
186
 
         if      ((e(i+1).gt.e(i)).and.(e(i).gt.e(i-1))) then
 
259
         if   ( ((e(i+1)-e(i)).gt.1.0e-2).and.
 
260
     >          ((e(i)-e(i-1)).gt.1.0d-2) ) then
187
261
            call dcopy(3*nion,t(rp),1,t(index),1)
188
 
         else if ((e(i-1).gt.e(i)).and.(e(i).gt.e(i+1))) then
 
262
         else if ( ((e(i-1)-e(i)).gt.1.0d-2) .and.
 
263
     >             ((e(i)-e(i+1)).gt.1.0d-2) ) then
189
264
            call dcopy(3*nion,t(rm),1,t(index),1)
190
265
         else
191
266
            
192
 
             if (dabs(e(i+1)-e(i)).gt.dabs(e(i-1)-e(i))) then
 
267
             if ( (dabs(e(i+1)-e(i))-dabs(e(i-1)-e(i)))
 
268
     >           .gt.1.0d-3) then
193
269
               dVmax = dabs(e(i+1)-e(i))
194
270
               dVmin = dabs(e(i-1)-e(i))
195
271
             else
197
273
               dVmin = dabs(e(i+1)-e(i))
198
274
             end if
199
275
 
200
 
             if (e(i+1).gt.e(i-1)) then
 
276
             if ((e(i+1)-e(i-1)).gt.1.0d-2) then
201
277
               call dscal(3*nion,dVmax,t(rp),1)
202
278
               call dscal(3*nion,dVmin,t(rm),1)
203
279
             else
222
298
      return
223
299
      end
224
300
 
 
301
*     ***********************************************
 
302
*     *                                             *
 
303
*     *             neb_gradient_get1               *
 
304
*     *                                             *
 
305
*     ***********************************************
225
306
 
226
307
*     *** RRR SUBRT TO GET AND ADJUST GRADIENT ***
227
308
*     *** RRR CHANGE TO INCLUDE CLIBING IMAGE ***
228
309
*     *** RRR CHANGE COMBINE PROJECTION and sforce
229
310
*     *** RRR addition into one subroutine
230
311
 
231
 
      subroutine neb_gradient_get1(bead_list,c,e,t,g)
 
312
      subroutine neb_gradient_get1(bead_list,kbeads,c,e,t,g)
232
313
 
233
314
 
234
315
      implicit none
235
316
      character*(*) bead_list
 
317
      real*8 kbeads
236
318
      real*8 c(*)
237
319
      real*8 e(*)
238
320
      real*8 t(*)
255
337
      external size_bead_list,nion_bead_list
256
338
      external ddot
257
339
 
258
 
*    *** RRR SET SPRING CONSTANT HERE
259
 
*    *** user-supplied k, default 0.1d0
260
 
*    *** currently not working
261
 
*    *** use default
262
 
*      ispring = rtdb_get(rtdb, 'spring', mt_dbl,1,spring)
263
 
*      write(*,*)'is spring constant set? =',ispring,spring
264
 
*      if(ispring) then
265
 
*         k = spring
266
 
*         write(*,*)'NEB spring constant set to  =',k
267
 
*      else
268
 
         k = 0.1d0
269
 
*         write(*,*)'NEB default springs =',k
270
 
*      endif      
271
 
*      stop
272
 
 
273
 
 
274
 
 
275
 
 
276
 
 
 
340
      k = kbeads
277
341
      nbeads = size_bead_list(bead_list)
278
342
      nion   = nion_bead_list(bead_list,1)
279
343
 
350
414
      return
351
415
      end
352
416
 
 
417
*     ***********************************************
 
418
*     *                                             *
 
419
*     *             neb_tangent1                    *
 
420
*     *                                             *
 
421
*     ***********************************************
353
422
 
354
423
*    *** RRR GET TANGENT TO BEAD PATH
355
424
*    *** RRR  DEFINE TANGENT AS THE VECTOR
431
500
      return
432
501
      end
433
502
 
 
503
*     ***********************************************
 
504
*     *                                             *
 
505
*     *             neb_energies_get                *
 
506
*     *                                             *
 
507
*     ***********************************************
 
508
 
434
509
      subroutine neb_energies_get(bead_list,e)
435
510
      implicit none
436
511
      character*(*) bead_list
454
529
      return
455
530
      end
456
531
 
 
532
*     ***********************************************
 
533
*     *                                             *
 
534
*     *             neb_coords_get                  *
 
535
*     *                                             *
 
536
*     ***********************************************
457
537
      subroutine neb_coords_get(bead_list,c)
458
538
      implicit none
459
539
      character*(*) bead_list
477
557
      return
478
558
      end
479
559
 
 
560
*     ***********************************************
 
561
*     *                                             *
 
562
*     *             neb_coords_set                  *
 
563
*     *                                             *
 
564
*     ***********************************************
480
565
      subroutine neb_coords_set(bead_list,c)
481
566
      implicit none
482
567
      character*(*) bead_list
496
581
        index = (i-1)*3*nion+1
497
582
        call coords_set_bead_list(bead_list,i,c(index))
498
583
      end do
499
 
 
500
 
 
501
584
      return
502
585
      end
503
586
 
 
587
*     ***********************************************
 
588
*     *                                             *
 
589
*     *             neb_masses_get                  *
 
590
*     *                                             *
 
591
*     ***********************************************
504
592
      subroutine  neb_masses_get(rtdb,m)
505
593
      implicit none
506
594
      integer rtdb
520
608
      external    bead_index_name
521
609
 
522
610
 
523
 
      geom_name   = 'neb_bead'//bead_index_name(1)//':geom'
 
611
      geom_name   = 'bead'//bead_index_name(1)//':geom'
524
612
      geomlen     = inp_strlen(geom_name)
525
613
 
526
614
      value = geom_create(geom,'neb_tmp')
532
620
      if (.not.value) call errquit('neb_masses failed',0,0)
533
621
 
534
622
      call dscal(nion,1822.89d0,m,1)
535
 
 
536
623
      return
537
624
      end
538
625
 
 
626
*     ***********************************************
 
627
*     *                                             *
 
628
*     *             neb_initialize                  *
 
629
*     *                                             *
 
630
*     ***********************************************
 
631
 
539
632
      subroutine neb_initialize(rtdb, bead_list)
540
633
      implicit none
541
634
      integer rtdb
552
645
c     This routine initializes the common /coptopt/ and
553
646
c     also creates and returns the geometry handle
554
647
c     
555
 
c      integer i, j, num, ma_type, nactive_atoms, l_actlist
556
648
      integer nbeads
557
 
c      logical ignore
558
 
      logical oprint
559
 
      character*80 title,neb_movecs
 
649
      character*80 neb_movecs
560
650
      logical custom_path
561
 
c     
562
 
      oprint = util_print('information', print_low)
563
 
     $     .and. (ga_nodeid() .eq. 0)
564
 
c
565
 
      if (oprint) then
566
 
         write(6,*)
567
 
         write(6,*)
568
 
         call util_print_centered(6,
569
 
     $        'NWChem Minimum Energy Pathway Program',
570
 
     $        40,.true.)
571
 
         write(6,*)
572
 
         write(6,*)
573
 
      endif
574
 
c
575
 
      if (rtdb_cget(rtdb,'title',1,title)) then
576
 
         if (oprint) then
577
 
            write(6,*)
578
 
            write(6,*)
579
 
            call util_print_centered(6, title, 40, .false.)
580
 
            write(6,*)
581
 
            write(6,*)
582
 
         endif
583
 
      endif
584
651
 
585
652
      if (.not.rtdb_cget(rtdb,'neb:movecs',1,neb_movecs)) then
586
653
         call util_file_prefix('movecs',neb_movecs)
592
659
      if (.not.rtdb_get(rtdb,'neb:nbeads',mt_int,1,nbeads)) then
593
660
         nbeads = 5
594
661
      end if
595
 
      if (oprint) then
596
 
         write(6,1) 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
597
 
     >              0,nbeads,
598
 
     >              neb_movecs
599
 
 1       format(
600
 
     $        ' maximum gradient threshold         (gmax) = ', f10.6,/,
601
 
     $        ' rms gradient threshold             (grms) = ', f10.6,/,
602
 
     $        ' maximum cartesian step threshold   (xmax) = ', f10.6,/,
603
 
     $        ' rms cartesian step threshold       (xrms) = ', f10.6,/,
604
 
     $        ' fixed trust radius                (trust) = ', f10.6,/,
605
 
     $        ' energy precision                  (eprec) = ', 1p,d9.1,
606
 
     $        0p,/,
607
 
     $        ' maximum number of steps          (nptopt) = ', i4,/,
608
 
     $        ' number of images in path         (nbeads) = ', i4,/,
609
 
     $        ' NEB movecs filename                       = ', a)
610
 
      end if
611
 
 
612
662
 
613
663
*     **** create bead_list *** 
614
664
      call init_bead_list(rtdb,bead_list,neb_movecs)
621
671
      return
622
672
      end
623
673
 
 
674
*     ***********************************************
 
675
*     *                                             *
 
676
*     *             neb_initial_path                *
 
677
*     *                                             *
 
678
*     ***********************************************
 
679
 
624
680
      subroutine neb_initial_path(rtdb,bead_list,nbeads)
625
681
      implicit none
626
682
      integer rtdb
627
683
      character*(*) bead_list
628
684
      integer nbeads
629
685
 
 
686
#include "geom.fh"
 
687
#include "rtdb.fh"
 
688
#include "util.fh"
 
689
#include "global.fh"
630
690
#include "mafdecls.fh"
631
 
#include "geom.fh"
 
691
#include "errquit.fh"
 
692
#include "stdio.fh"
632
693
 
633
694
*     **** local variables ****
634
 
      logical value
635
 
      integer i,geom,geomlen,movecslen,nion
636
 
      integer r1(2),r2(2),r3(2)
637
 
      real*8  t
 
695
      logical value,impose,oprint,hasmiddle
 
696
      integer i,geom,geomlen,movecslen,nion,nfit,i2,i3,pathguess
 
697
      integer r1(2),r2(2),r3(2),rmid(2),ifit(2),wfit(2),rcoords(2)
 
698
      real*8  t,rms1,rms2
638
699
      character*255 geom_name,movecs_name
639
700
 
640
701
*     **** external functions ****
643
704
      external    inp_strlen
644
705
      external    bead_index_name
645
706
 
 
707
      oprint = (ga_nodeid().eq.0)
 
708
 
 
709
      if (.not. rtdb_get(rtdb,'neb:pathguess',mt_int,1,pathguess))
 
710
     $      pathguess = 2
 
711
 
646
712
      value = geom_create(geom,'neb_tmp')
647
 
      value = value.and.geom_rtdb_load(rtdb,geom,'neb_start')
 
713
 
 
714
      !*** check for neb_start, otherwise just read geometry ****
 
715
      if (.not.geom_rtdb_load(rtdb,geom,'neb_start')) then
 
716
         value = value.and.geom_rtdb_load(rtdb,geom,'geometry')
 
717
      end if
 
718
 
648
719
      value = value.and.geom_ncent(geom,nion)
649
720
      value = value.and.geom_destroy(geom)
650
721
      if (.not.value) call errquit('neb_initial_path failed',0,0)
652
723
      value = value.and.MA_push_get(mt_dbl, (3*nion), 'r1',r1(2),r1(1))
653
724
      value = value.and.MA_push_get(mt_dbl, (3*nion), 'r2',r2(2),r2(1))
654
725
      value = value.and.MA_push_get(mt_dbl, (3*nion), 'r3',r3(2),r3(1))
 
726
      value = value.and.MA_push_get(mt_dbl, (3*nion),
 
727
     >                              'rmid',rmid(2),rmid(1))
655
728
      if (.not.value) call errquit('neb_initial_path failed',1,0)
656
729
 
657
730
      value = value.and.geom_create(geom,'neb_tmp')
658
 
      value = value.and.geom_rtdb_load(rtdb,geom,'neb_end')
 
731
 
 
732
      !*** try neb_end, then endgeom ****
 
733
      if (.not.geom_rtdb_load(rtdb,geom,'neb_end')) then
 
734
         value = value.and.geom_rtdb_load(rtdb,geom,'endgeom')
 
735
      end if
659
736
      value = value.and.geom_cart_coords_get(geom,dbl_mb(r2(1)))
660
 
      value = value.and.geom_rtdb_load(rtdb,geom,'neb_start')
 
737
 
 
738
      !*** try neb_start, then geometry****
 
739
      if (.not.geom_rtdb_load(rtdb,geom,'neb_start')) then
 
740
         value = value.and.geom_rtdb_load(rtdb,geom,'geometry')
 
741
      end if
661
742
      value = value.and.geom_cart_coords_get(geom,dbl_mb(r1(1)))
662
743
      if (.not.value) call errquit('neb_initial_path failed',2,0)
663
744
 
664
 
      do i=1,nbeads
665
 
        t = (i-1)/dble(nbeads-1)
666
 
 
667
 
        call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(r3(1)),1)
668
 
        call dscal(3*nion,(1.0d0-t),dbl_mb(r3(1)),1)
669
 
        call daxpy(3*nion,t,dbl_mb(r2(1)),1,dbl_mb(r3(1)),1)
670
 
 
671
 
        geom_name   = 'neb_bead'//bead_index_name(i)//':geom'
672
 
        movecs_name = 'neb_bead'//bead_index_name(i)//'.movecs'
673
 
        geomlen     = inp_strlen(geom_name)
674
 
        movecslen   = inp_strlen(movecs_name)
675
 
        value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
676
 
        value = value.and.geom_rtdb_store(rtdb,geom,
677
 
     >                                    geom_name(1:geomlen))
678
 
 
679
 
        call add_bead_list(bead_list,
680
 
     >                     movecs_name(1:movecslen),
681
 
     >                     geom_name(1:geomlen))
682
 
 
683
 
      end do
 
745
      if (.not.rtdb_get(rtdb,'neb:hasmiddle',mt_log,1,hasmiddle))
 
746
     >   hasmiddle = .false.
 
747
 
 
748
      if (hasmiddle) then
 
749
         if (.not.geom_rtdb_load(rtdb,geom,'neb_middle')) then
 
750
            if (.not.geom_rtdb_load(rtdb,geom,'midgeom')) then
 
751
               hasmiddle = .false.
 
752
            end if
 
753
         else
 
754
            if (.not.geom_cart_coords_get(geom,dbl_mb(rmid(1))))
 
755
     >      call errquit('neb_initial_path failed',2,0)
 
756
         end if
 
757
      end if
 
758
 
 
759
      if (oprint) then
 
760
         write(luout,*) 
 
761
     >   " - Generating initial path by linear interpolation"
 
762
         write(luout,'(A,I4)') "    + number images = ",nbeads
 
763
         if (hasmiddle) then
 
764
            write(luout,'(3A)')    
 
765
     >      "    + neb_start (geometry) geometry -->",
 
766
     >      " neb_middle (midgeom) geometry -->",
 
767
     >      " neb_end (endgeom) geometry"
 
768
         else
 
769
            write(luout,'(A)')    
 
770
     >      "    + neb_start (geometry) geometry -->",
 
771
     >      " neb_end (endgeom) geometry"
 
772
         end if
 
773
      end if
 
774
 
 
775
 
 
776
 
 
777
      if (.not.rtdb_get(rtdb,'neb:impose',mt_log,1,impose))
 
778
     >   impose = .false.
 
779
 
 
780
      if (impose) then
 
781
         value =  MA_push_get(mt_int,(2*nion),'ifit',ifit(2),ifit(1))
 
782
         value = value.and.
 
783
     >            MA_push_get(mt_dbl,(nion),'wfit',wfit(2),wfit(1))
 
784
         if(.not.value) call errquit('neb_initial_path failed',3,MA_ERR)
 
785
 
 
786
         if (hasmiddle) then
 
787
            call neb_impose(nion,dbl_mb(r1(1)),dbl_mb(rmid(1)),
 
788
     >                   nfit,int_mb(ifit(1)),dbl_mb(wfit(1)),rms1,rms2)
 
789
            if (oprint) then
 
790
            write(luout,*) 
 
791
     >      " - Imposing neb_mid (midgeom) geometry onto",
 
792
     >      " neb_start (geometry) geometry"
 
793
            write(luout,'(A,F10.6)') "    + initial rmsq = ",rms1
 
794
            write(luout,'(A,F10.6)') "    + imposed rmsq = ",rms2
 
795
            end if
 
796
 
 
797
            call neb_impose(nion,dbl_mb(rmid(1)),dbl_mb(r2(1)),
 
798
     >                   nfit,int_mb(ifit(1)),dbl_mb(wfit(1)),rms1,rms2)
 
799
            if (oprint) then
 
800
            write(luout,*) 
 
801
     >      " - Imposing neb_end geometry onto neb_mid geometry"
 
802
            write(luout,'(A,F10.6)') "    + initial rmsq = ",rms1
 
803
            write(luout,'(A,F10.6)') "    + imposed rmsq = ",rms2
 
804
            end if
 
805
         else
 
806
            call neb_impose(nion,dbl_mb(r1(1)),dbl_mb(r2(1)),
 
807
     >                   nfit,int_mb(ifit(1)),dbl_mb(wfit(1)),rms1,rms2)
 
808
            if (oprint) then
 
809
            write(luout,*) 
 
810
     >      " - Imposing neb_end (endgeom) geometry",
 
811
     >      " onto neb_start (geometry) geometry"
 
812
            write(luout,'(A,F10.6)') "    + initial rmsq = ",rms1
 
813
            write(luout,'(A,F10.6)') "    + imposed rmsq = ",rms2
 
814
            end if
 
815
         end if
 
816
 
 
817
         value =           MA_pop_stack(wfit(2))
 
818
         value = value.and.MA_pop_stack(ifit(2))
 
819
         if(.not.value) call errquit('neb_initial_path failed',4,MA_ERR)
 
820
      end if
 
821
 
 
822
 
 
823
 
 
824
      if(.not.ma_push_get(mt_dbl,3*nion*nbeads,'rcoords',
 
825
     >                    rcoords(2), rcoords(1)))
 
826
     >     call errquit('neb_initialize: memory',3*nion*nbeads,MA_ERR)
 
827
 
 
828
 
 
829
 
 
830
      if (hasmiddle) then
 
831
         i2 = nbeads/2+1
 
832
 
 
833
        call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(rcoords(1)),1) 
 
834
        call dcopy(3*nion,dbl_mb(rmid(1)),1,
 
835
     >            dbl_mb(rcoords(1)+3*nion*(i2-1)),1) 
 
836
        call dcopy(3*nion,dbl_mb(r2(1)),1,
 
837
     >            dbl_mb(rcoords(1)+3*nion*(nbeads-1)),1) 
 
838
 
 
839
        call zts_guess(nion,nbeads,i2,dbl_mb(rcoords(1)),pathguess)
 
840
 
 
841
        do i=1,nbeads
 
842
          call dcopy(3*nion,dbl_mb(rcoords(1)+3*nion*(i-1)),1,
 
843
     >                      dbl_mb(r3(1)),1)
 
844
 
 
845
          geom_name   = 'bead'//bead_index_name(i)//':geom'
 
846
          movecs_name = 'bead'//bead_index_name(i)//'.movecs'
 
847
          geomlen     = inp_strlen(geom_name)
 
848
          movecslen   = inp_strlen(movecs_name)
 
849
          value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
 
850
          value = value.and.geom_rtdb_store(rtdb,geom,
 
851
     >                                    geom_name(1:geomlen))
 
852
          call add_bead_list(bead_list,
 
853
     >                     movecs_name(1:movecslen),
 
854
     >                     geom_name(1:geomlen))
 
855
        end do
 
856
 
 
857
 
 
858
c        do i=1,i2
 
859
c          t = (i-1)/dble(i2-1)
 
860
c          call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(r3(1)),1)
 
861
c          call dscal(3*nion,(1.0d0-t),dbl_mb(r3(1)),1)
 
862
c          call daxpy(3*nion,t,dbl_mb(rmid(1)),1,dbl_mb(r3(1)),1)
 
863
c
 
864
c          geom_name   = 'bead'//bead_index_name(i)//':geom'
 
865
c          movecs_name = 'bead'//bead_index_name(i)//'.movecs'
 
866
c          geomlen     = inp_strlen(geom_name)
 
867
c          movecslen   = inp_strlen(movecs_name)
 
868
c          value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
 
869
c          value = value.and.geom_rtdb_store(rtdb,geom,
 
870
c     >                                    geom_name(1:geomlen))
 
871
c          call add_bead_list(bead_list,
 
872
c     >                     movecs_name(1:movecslen),
 
873
c     >                     geom_name(1:geomlen))
 
874
c        end do
 
875
c
 
876
c        i3 = nbeads-i2
 
877
c        do i=1,i3
 
878
c          t = i/dble(i3)
 
879
c          call dcopy(3*nion,dbl_mb(rmid(1)),1,dbl_mb(r3(1)),1)
 
880
c          call dscal(3*nion,(1.0d0-t),dbl_mb(r3(1)),1)
 
881
c          call daxpy(3*nion,t,dbl_mb(r2(1)),1,dbl_mb(r3(1)),1)
 
882
c
 
883
c          geom_name   = 'bead'//bead_index_name(i+i2)//':geom'
 
884
c          movecs_name = 'bead'//bead_index_name(i+i2)//'.movecs'
 
885
c          geomlen     = inp_strlen(geom_name)
 
886
c          movecslen   = inp_strlen(movecs_name)
 
887
c          value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
 
888
c          value = value.and.geom_rtdb_store(rtdb,geom,
 
889
c     >                                    geom_name(1:geomlen))
 
890
c          call add_bead_list(bead_list,
 
891
c     >                     movecs_name(1:movecslen),
 
892
c     >                     geom_name(1:geomlen))
 
893
c        end do
 
894
 
 
895
*     **** linear interpolation with Robinson Checking ****
 
896
      else
 
897
      
 
898
        call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(rcoords(1)),1) 
 
899
        call dcopy(3*nion,dbl_mb(r2(1)),1,
 
900
     >            dbl_mb(rcoords(1)+3*nion*(nbeads-1)),1) 
 
901
 
 
902
        call zts_guessall(nion,nbeads,dbl_mb(rcoords(1)),geom)
 
903
 
 
904
        do i=1,nbeads
 
905
c          t = (i-1)/dble(nbeads-1)
 
906
c
 
907
c          call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(r3(1)),1)
 
908
c          call dscal(3*nion,(1.0d0-t),dbl_mb(r3(1)),1)
 
909
c          call daxpy(3*nion,t,dbl_mb(r2(1)),1,dbl_mb(r3(1)),1)
 
910
 
 
911
          call dcopy(3*nion,dbl_mb(rcoords(1)+3*nion*(i-1)),1,
 
912
     >                      dbl_mb(r3(1)),1)
 
913
 
 
914
          geom_name   = 'bead'//bead_index_name(i)//':geom'
 
915
          movecs_name = 'bead'//bead_index_name(i)//'.movecs'
 
916
          geomlen     = inp_strlen(geom_name)
 
917
          movecslen   = inp_strlen(movecs_name)
 
918
          value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
 
919
          value = value.and.geom_rtdb_store(rtdb,geom,
 
920
     >                                    geom_name(1:geomlen))
 
921
 
 
922
          call add_bead_list(bead_list,
 
923
     >                     movecs_name(1:movecslen),
 
924
     >                     geom_name(1:geomlen))
 
925
 
 
926
        end do
 
927
      end if
 
928
      value = value.and.MA_pop_stack(rcoords(2))
684
929
      value = value.and.geom_destroy(geom)
 
930
      value = value.and.MA_pop_stack(rmid(2))
685
931
      value = value.and.MA_pop_stack(r3(2))
686
932
      value = value.and.MA_pop_stack(r2(2))
687
933
      value = value.and.MA_pop_stack(r1(2))
691
937
      return
692
938
      end
693
939
 
 
940
*     ***********************************************
 
941
*     *                                             *
 
942
*     *             neb_initial_path_custom         *
 
943
*     *                                             *
 
944
*     ***********************************************
694
945
      subroutine neb_initial_path_custom(rtdb,bead_list,nbeads)
695
946
      implicit none
696
947
      integer rtdb
712
963
      external    bead_index_name
713
964
 
714
965
      do i=1,nbeads
715
 
        geom_name   = 'neb_bead'//bead_index_name(i)//':geom'
716
 
        movecs_name = 'neb_bead'//bead_index_name(i)//'.movecs'
 
966
        geom_name   = 'bead'//bead_index_name(i)//':geom'
 
967
        movecs_name = 'bead'//bead_index_name(i)//'.movecs'
717
968
        geomlen     = inp_strlen(geom_name)
718
969
        movecslen   = inp_strlen(movecs_name)
719
970
 
726
977
      return
727
978
      end
728
979
 
 
980
*     ***********************************************
 
981
*     *                                             *
 
982
*     *             neb_verlet_update               *
 
983
*     *                                             *
 
984
*     ***********************************************
 
985
 
729
986
      subroutine neb_verlet_update(ng,c0,c1,v1,dti,g1)
730
987
        integer ng
731
988
        double precision c0(*)
756
1013
      end
757
1014
 
758
1015
      subroutine neb_cg_direction(ng,g0,g1,s)
759
 
        integer ng
760
 
        double precision g0(*)
761
 
        double precision g1(*)
762
 
        double precision s(*)
763
 
 
764
 
        integer i
765
 
        double precision gamma1
766
 
        double precision gamma2
767
 
        double precision sn
768
 
 
769
 
c       *** choosing Polac-Ribiere coeff ***        
770
 
        gamma1=0.0d0
771
 
        gamma2=0.0d0
772
 
        do i=1,ng
773
 
           gamma1 = gamma1 + (g1(i)-g0(i))*g1(i)
774
 
           gamma2 = gamma2 + g0(i)*g0(i)
775
 
        end do
776
 
 
777
 
        do i=1,ng
778
 
           s(i) = -g1(i) + s(i)*gamma1/gamma2
779
 
        end do
 
1016
      integer ng
 
1017
      double precision g0(*)
 
1018
      double precision g1(*)
 
1019
      double precision s(*)
 
1020
 
 
1021
      integer i
 
1022
      double precision gamma1
 
1023
      double precision gamma2
 
1024
      double precision sn
 
1025
 
 
1026
c     *** choosing Polac-Ribiere coeff ***        
 
1027
      gamma1=0.0d0
 
1028
      gamma2=0.0d0
 
1029
      do i=1,ng
 
1030
         gamma1 = gamma1 + (g1(i)-g0(i))*g1(i)
 
1031
         gamma2 = gamma2 + g0(i)*g0(i)
 
1032
      end do
 
1033
 
 
1034
      do i=1,ng
 
1035
         s(i) = -g1(i) + s(i)*gamma1/gamma2
 
1036
      end do
780
1037
        
 
1038
      return
781
1039
      end
782
1040
 
783
1041
      subroutine neb_move(ng,dt,c0,c1,s)
784
 
        integer ng
785
 
        double precision dt
786
 
        double precision c0(*)
787
 
        double precision c1(*)
788
 
        double precision s(*)
789
 
 
790
 
        integer i
791
 
 
792
 
 
793
 
        do i=1,ng
794
 
            c1(i) = c0(i)+dt*s(i)
795
 
        end do
796
 
 
 
1042
      integer ng
 
1043
      double precision dt
 
1044
      double precision c0(*)
 
1045
      double precision c1(*)
 
1046
      double precision s(*)
 
1047
 
 
1048
      integer i
 
1049
 
 
1050
      do i=1,ng
 
1051
          c1(i) = c0(i)+dt*s(i)
 
1052
      end do
 
1053
 
 
1054
      return
797
1055
      end
798
1056
 
799
1057
      subroutine neb_calc_convergence(ng,g1,c0,c1,Gmax,Grms,Xmax,Xrms)
804
1062
      double precision c1(*)
805
1063
      double precision Gmax,Grms,Xmax,Xrms
806
1064
 
 
1065
#include "global.fh"
 
1066
 
807
1067
      !***** local variables ****
808
 
      integer ii
 
1068
      logical oprint
 
1069
      integer ii,imax
809
1070
      double precision dx
810
1071
 
 
1072
 
 
1073
      oprint = ga_nodeid() .eq. 0
 
1074
 
811
1075
      Gmax = 0.0
 
1076
      imax = 0
812
1077
      do ii=1,ng
813
 
        if (dabs(g1(ii)).gt.Gmax) Gmax = dabs(g1(ii))
 
1078
        if (dabs(g1(ii)).gt.Gmax) then
 
1079
            Gmax = dabs(g1(ii))
 
1080
            imax = ii
 
1081
        end if
814
1082
      end do
 
1083
      if (oprint) write(*,*) "neb: imax,Gmax=",imax,Gmax
815
1084
 
816
1085
      Grms = 0.0
817
1086
      do ii=1,ng
835
1104
 
836
1105
      return
837
1106
      end
 
1107
 
 
1108
 
 
1109
*     *****************************************
 
1110
*     *                                       *
 
1111
*     *             neb_lmbfgs                *
 
1112
*     *                                       *
 
1113
*     *****************************************
 
1114
      subroutine neb_lmbfgs(n,m,x,g,hg)
 
1115
      implicit none
 
1116
      integer n,m
 
1117
      real*8 x(n,m)
 
1118
      real*8 g(n,m)
 
1119
      real*8 hg(n)
 
1120
 
 
1121
*     **** local variables ****
 
1122
      integer k
 
1123
      real*8 rho(25),alpha(25),beta(25)
 
1124
      real*8 tmp
 
1125
 
 
1126
*     **** external functions ****
 
1127
      real*8   ddot
 
1128
      external ddot
 
1129
 
 
1130
*     **** compute rho(k) = 1/y(:,k)' * s(:,k) ****
 
1131
      do k=1,m-1
 
1132
         tmp =       ddot(n,x(1,k+1),1,g(1,k+1),1)
 
1133
         tmp = tmp - ddot(n,x(1,k+1),1,g(1,k),  1)
 
1134
         tmp = tmp - ddot(n,x(1,k),  1,g(1,k+1),1)
 
1135
         tmp = tmp + ddot(n,x(1,k),  1,g(1,k),  1)
 
1136
         if (dabs(tmp).gt.1.0d-9) then
 
1137
             rho(k) = 1.0d0/tmp
 
1138
         else
 
1139
             rho(k) = 0.0d0
 
1140
         end if
 
1141
      end do
 
1142
 
 
1143
      call dcopy(n,g(1,m),1,hg,1)
 
1144
 
 
1145
      do k = (m-1),1,-1
 
1146
        alpha(k) = rho(k)
 
1147
     >            *(ddot(n,x(1,k+1),1,hg,1) - ddot(n,x(1,k),1,hg,1))
 
1148
 
 
1149
        call daxpy(n,(-alpha(k)),g(1,k+1),1,hg,1)
 
1150
        call daxpy(n,( alpha(k)),g(1,k),  1,hg,1)
 
1151
      end do
 
1152
 
 
1153
      do k = 1,(m-1)
 
1154
        beta(k) = rho(k)
 
1155
     >           *(ddot(n,g(1,k+1),1,hg,1) - ddot(n,g(1,k),1,hg,1))
 
1156
 
 
1157
        call daxpy(n,(alpha(k)-beta(k)),x(1,k+1),1,hg,1)
 
1158
        call daxpy(n,(beta(k)-alpha(k)),x(1,k),  1,hg,1)
 
1159
      end do
 
1160
 
 
1161
      return
 
1162
      end
 
1163
 
 
1164
 
 
1165
*     ***********************************************
 
1166
*     *                                             *
 
1167
*     *             neb_resize_path                 *
 
1168
*     *                                             *
 
1169
*     ***********************************************
 
1170
      subroutine neb_resize_path(rtdb,bead_list,nbeads1,nbeads2)
 
1171
      implicit none
 
1172
      integer rtdb
 
1173
      character*(*) bead_list
 
1174
      integer nbeads1
 
1175
      integer nbeads2
 
1176
 
 
1177
#include "mafdecls.fh"
 
1178
#include "geom.fh"
 
1179
 
 
1180
*     **** local variables ****
 
1181
      logical value
 
1182
      integer i,geom,geomlen,movecslen,nion
 
1183
      integer j1,j2,shift
 
1184
      integer r1(2),r2(2),r3(2),c(2)
 
1185
      real*8  t,t1,t2,t3
 
1186
      character*255 geom_name,movecs_name
 
1187
 
 
1188
*     **** external functions ****
 
1189
      integer     inp_strlen
 
1190
      character*7 bead_index_name
 
1191
      external    inp_strlen
 
1192
      external    bead_index_name
 
1193
 
 
1194
      value = geom_create(geom,'neb_tmp')
 
1195
      if (.not.geom_rtdb_load(rtdb,geom,'neb_start')) then
 
1196
          value = value.and.geom_rtdb_load(rtdb,geom,'geometry')
 
1197
      end if
 
1198
      value = value.and.geom_ncent(geom,nion)
 
1199
      value = value.and.geom_destroy(geom)
 
1200
      if (.not.value) call errquit('neb_resize_path failed',0,0)
 
1201
 
 
1202
 
 
1203
      value = value.and.MA_push_get(mt_dbl,(3*nion*nbeads1),
 
1204
     >                              'c',c(2),c(1))
 
1205
      value = value.and.MA_push_get(mt_dbl,(3*nion),'r1',r1(2),r1(1))
 
1206
      value = value.and.MA_push_get(mt_dbl,(3*nion),'r2',r2(2),r2(1))
 
1207
      value = value.and.MA_push_get(mt_dbl,(3*nion),'r3',r3(2),r3(1))
 
1208
      if (.not.value) call errquit('neb_resize_path failed',1,0)
 
1209
 
 
1210
      value = value.and.geom_create(geom,'neb_tmp')
 
1211
      if (.not.geom_rtdb_load(rtdb,geom,'neb_start')) then
 
1212
         value = value.and.geom_rtdb_load(rtdb,geom,'geometry')
 
1213
      end if
 
1214
      value = value.and.geom_cart_coords_get(geom,dbl_mb(r3(1)))
 
1215
      if (.not.value) call errquit('neb_resize_path failed',2,0)
 
1216
 
 
1217
      do i=1,nbeads1
 
1218
         shift = (i-1)*3*nion
 
1219
         call coords_get_bead_list(bead_list,i,dbl_mb(c(1)+shift))
 
1220
      end do
 
1221
      call reset_bead_list(bead_list)
 
1222
 
 
1223
      do i=1,nbeads2
 
1224
        t = (i-1)/dble(nbeads2-1)
 
1225
 
 
1226
        j1 = t*(nbeads1-1) + 1
 
1227
        j2 = j1+1
 
1228
        t1 = (j1-1)/dble(nbeads1-1)
 
1229
        t2 = (j2-1)/dble(nbeads1-1)
 
1230
        t3 = (t-t1)/(t2-t1)
 
1231
 
 
1232
        if (j2.gt.nbeads1) then
 
1233
           t3 = 0.0d0
 
1234
           j2=nbeads1
 
1235
        end if
 
1236
 
 
1237
        shift = (j1-1)*3*nion
 
1238
        call dcopy(3*nion,dbl_mb(c(1)+shift),1,dbl_mb(r1(1)),1)
 
1239
 
 
1240
        shift = (j2-1)*3*nion
 
1241
        call dcopy(3*nion,dbl_mb(c(1)+shift),1,dbl_mb(r2(1)),1)
 
1242
 
 
1243
        call dcopy(3*nion,dbl_mb(r1(1)),1,dbl_mb(r3(1)),1)
 
1244
        call dscal(3*nion,(1.0d0-t3),dbl_mb(r3(1)),1)
 
1245
        call daxpy(3*nion,t3,dbl_mb(r2(1)),1,dbl_mb(r3(1)),1)
 
1246
 
 
1247
        geom_name   = 'bead'//bead_index_name(i)//':geom'
 
1248
        movecs_name = 'bead'//bead_index_name(i)//'.movecs'
 
1249
        geomlen     = inp_strlen(geom_name)
 
1250
        movecslen   = inp_strlen(movecs_name)
 
1251
        value = value.and.geom_cart_coords_set(geom,dbl_mb(r3(1)))
 
1252
        value = value.and.geom_rtdb_store(rtdb,geom,
 
1253
     >                                    geom_name(1:geomlen))
 
1254
 
 
1255
        call add_bead_list(bead_list,
 
1256
     >                     movecs_name(1:movecslen),
 
1257
     >                     geom_name(1:geomlen))
 
1258
 
 
1259
      end do
 
1260
      value = value.and.geom_destroy(geom)
 
1261
      value = value.and.MA_pop_stack(r3(2))
 
1262
      value = value.and.MA_pop_stack(r2(2))
 
1263
      value = value.and.MA_pop_stack(r1(2))
 
1264
      value = value.and.MA_pop_stack(c(2))
 
1265
      if (.not.value) call errquit('neb_new_path failed',3,0)
 
1266
      return
 
1267
      end
 
1268
 
 
1269
 
 
1270
*     ***********************************************
 
1271
*     *                                             *
 
1272
*     *             neb_impose                      *
 
1273
*     *                                             *
 
1274
*     ***********************************************
 
1275
c     subroutine neb_impose  --  superimpose two coordinate sets
 
1276
c
 
1277
c     This routine performs the least squares best superposition
 
1278
c     of two atomic coordinate sets via a quaternion method;
 
1279
c     upon return, the first coordinate set is unchanged while
 
1280
c     the second set is translated and rotated to give best fit;
 
1281
c     the final root mean square fit is returned in "rmsvalue"
 
1282
c
 
1283
      subroutine neb_impose(nion,rion1,rion2,nfit,ifit,wfit,rms1,rms2)
 
1284
      implicit none
 
1285
      integer nion
 
1286
      real*8 rion1(3,*),rion2(3,*)
 
1287
      integer nfit,ifit(2,*)
 
1288
      real*8  wfit(*)
 
1289
      real*8 rms1,rms2
 
1290
 
 
1291
*     **** local variables ****
 
1292
      integer i
 
1293
      real*8 xmid,ymid,zmid
 
1294
 
 
1295
*     **** external functions ****
 
1296
      real*8   neb_rmsfit
 
1297
      external neb_rmsfit
 
1298
 
 
1299
      nfit = nion
 
1300
      do i=1,nfit
 
1301
         ifit(1,i) = i
 
1302
         ifit(2,i) = i
 
1303
         wfit(i)   = 1.0d0
 
1304
      end do
 
1305
      rms1 = neb_rmsfit(nfit,ifit,wfit,rion1,rion2)
 
1306
c
 
1307
c     superimpose the centroids of active atom pairs
 
1308
c
 
1309
      call neb_center(nion,rion1,rion2,
 
1310
     >                      nfit,ifit,wfit,
 
1311
     >                      xmid,ymid,zmid)
 
1312
 
 
1313
c
 
1314
c     use a quaternion method to achieve the superposition
 
1315
c
 
1316
      call neb_quatfit(nion,rion1,rion2,nfit,ifit,wfit)
 
1317
c
 
1318
c     translate both coordinate sets so as to return
 
1319
c     the first set to its original position
 
1320
c
 
1321
      do i = 1, nion
 
1322
         rion1(1,i) = rion1(1,i) + xmid
 
1323
         rion1(2,i) = rion1(2,i) + ymid
 
1324
         rion1(3,i) = rion1(3,i) + zmid
 
1325
      end do
 
1326
      do i = 1, nion
 
1327
         rion2(1,i) = rion2(1,i) + xmid
 
1328
         rion2(2,i) = rion2(2,i) + ymid
 
1329
         rion2(3,i) = rion2(3,i) + zmid
 
1330
      end do
 
1331
      rms2 = neb_rmsfit(nfit,ifit,wfit,rion1,rion2)
 
1332
 
 
1333
      return
 
1334
      end
 
1335
 
 
1336
*     ***********************************************
 
1337
*     *                                             *
 
1338
*     *             neb_rmsfit                      *
 
1339
*     *                                             *
 
1340
*     ***********************************************
 
1341
c     function neb_rmsfit  --  rms deviation for paired atoms
 
1342
c
 
1343
c     This routine computes the rms fit of two coordinate sets
 
1344
c
 
1345
      real*8 function neb_rmsfit(nfit,ifit,wfit,rion1,rion2)
 
1346
      implicit none
 
1347
      integer nfit,ifit(2,*)
 
1348
      real*8 wfit(*)
 
1349
      real*8 rion1(3,*),rion2(3,*)
 
1350
 
 
1351
      integer i,i1,i2
 
1352
      real*8 rmsterm,rmsfit
 
1353
      real*8 xr,yr,zr,dist2
 
1354
      real*8 weight,norm
 
1355
c
 
1356
c     compute the rms fit over superimposed atom pairs
 
1357
c
 
1358
      rmsfit = 0.0d0
 
1359
      norm = 0.0d0
 
1360
      do i = 1, nfit
 
1361
         i1 = ifit(1,i)
 
1362
         i2 = ifit(2,i)
 
1363
         weight = wfit(i)
 
1364
         xr = rion1(1,i1) - rion2(1,i2)
 
1365
         yr = rion1(2,i1) - rion2(2,i2)
 
1366
         zr = rion1(3,i1) - rion2(3,i2)
 
1367
         dist2 = xr**2 + yr**2 + zr**2
 
1368
         norm = norm + weight
 
1369
         rmsterm = dist2 * weight
 
1370
         rmsfit = rmsfit + rmsterm
 
1371
      end do
 
1372
      neb_rmsfit = sqrt(rmsfit/norm)
 
1373
      return
 
1374
      end
 
1375
 
 
1376
 
 
1377
*     ***********************************************
 
1378
*     *                                             *
 
1379
*     *             neb_center                      *
 
1380
*     *                                             *
 
1381
*     ***********************************************
 
1382
c     subroutine neb_center  --  superimpose structure centroids 
 
1383
c
 
1384
c     This routine moves the weighted centroid of each coordinate
 
1385
c     set to the origin during least squares superposition
 
1386
 
 
1387
      subroutine neb_center(nion,rion1,rion2,
 
1388
     >                      nfit,ifit,wfit,
 
1389
     >                      xmid,ymid,zmid)
 
1390
      implicit none
 
1391
      integer nion
 
1392
      real*8 rion1(3,*),rion2(3,*)
 
1393
      integer nfit,ifit(2,*)
 
1394
      real*8  wfit(*)
 
1395
      real*8 xmid,ymid,zmid
 
1396
 
 
1397
*     **** local variables ****
 
1398
      integer i,k
 
1399
      real*8 weight,norm
 
1400
c
 
1401
c
 
1402
c     find the weighted centroid of the second
 
1403
c     structure and translate it to the origin
 
1404
c
 
1405
      xmid = 0.0d0
 
1406
      ymid = 0.0d0
 
1407
      zmid = 0.0d0
 
1408
      norm = 0.0d0
 
1409
      do i = 1, nfit
 
1410
         k = ifit(2,i)
 
1411
         weight = wfit(i)
 
1412
         xmid = xmid + rion2(1,k)*weight
 
1413
         ymid = ymid + rion2(2,k)*weight
 
1414
         zmid = zmid + rion2(3,k)*weight
 
1415
         norm = norm + weight
 
1416
      end do
 
1417
 
 
1418
      xmid = xmid / norm
 
1419
      ymid = ymid / norm
 
1420
      zmid = zmid / norm
 
1421
      do i = 1, nion
 
1422
         rion2(1,i) = rion2(1,i) - xmid
 
1423
         rion2(2,i) = rion2(2,i) - ymid
 
1424
         rion2(3,i) = rion2(3,i) - zmid
 
1425
      end do
 
1426
c
 
1427
c     now repeat for the first structure, note
 
1428
c     that this centroid position gets returned
 
1429
c
 
1430
      xmid = 0.0d0
 
1431
      ymid = 0.0d0
 
1432
      zmid = 0.0d0
 
1433
      norm = 0.0d0
 
1434
      do i = 1, nfit
 
1435
         k = ifit(1,i)
 
1436
         weight = wfit(i)
 
1437
         xmid = xmid + rion1(1,k)*weight
 
1438
         ymid = ymid + rion1(2,k)*weight
 
1439
         zmid = zmid + rion1(3,k)*weight
 
1440
         norm = norm + weight
 
1441
      end do
 
1442
 
 
1443
      xmid = xmid / norm
 
1444
      ymid = ymid / norm
 
1445
      zmid = zmid / norm
 
1446
      do i = 1, nion
 
1447
         rion1(1,i) = rion1(1,i) - xmid
 
1448
         rion1(2,i) = rion1(2,i) - ymid
 
1449
         rion1(3,i) = rion1(3,i) - zmid
 
1450
      end do
 
1451
 
 
1452
      return
 
1453
      end
 
1454
 
 
1455
 
 
1456
*     ***********************************************
 
1457
*     *                                             *
 
1458
*     *             neb_quatfit                     *
 
1459
*     *                                             *
 
1460
*     ***********************************************
 
1461
c     subroutine quatfit  --  quaternion superposition of coords
 
1462
c
 
1463
c     This routine uses a quaternion-based method to achieve the best
 
1464
c     fit superposition of two sets of coordinates
 
1465
c
 
1466
c     literature reference:
 
1467
c
 
1468
c     S. J. Kearsley, "An Algorithm for the Simultaneous Superposition
 
1469
c     of a Structural Series", Journal of Computational Chemistry,
 
1470
c     11, 1187-1192 (1990)
 
1471
c
 
1472
c     adapted from an original program written by David J. Heisterberg,
 
1473
c     Ohio Supercomputer Center, Columbus, OH
 
1474
c
 
1475
      subroutine neb_quatfit(nion,rion1,rion2,nfit,ifit,wfit)
 
1476
      implicit none
 
1477
      integer nion
 
1478
      real*8 rion1(3,*),rion2(3,*)
 
1479
      integer nfit
 
1480
      integer ifit(2,*)
 
1481
      real*8 wfit(*)
 
1482
 
 
1483
 
 
1484
      integer i,i1,i2,n1,n2
 
1485
      real*8  weight,xrot,drot,zrot
 
1486
      real*8  xxyx,xxyy,xxyz,xyyx,xyyy
 
1487
      real*8  xyyz,xzyx,xzyy,xzyz
 
1488
      real*8  rot(3,3),temp1(4),temp2(4)
 
1489
      real*8  q(4),d(4),c(4,4),v(4,4)
 
1490
c
 
1491
c     build the upper triangle of the quadratic form matrix
 
1492
c
 
1493
      xxyx = 0.0d0
 
1494
      xxyy = 0.0d0
 
1495
      xxyz = 0.0d0
 
1496
      xyyx = 0.0d0
 
1497
      xyyy = 0.0d0
 
1498
      xyyz = 0.0d0
 
1499
      xzyx = 0.0d0
 
1500
      xzyy = 0.0d0
 
1501
      xzyz = 0.0d0
 
1502
      do i = 1, nfit
 
1503
         i1 = ifit(1,i)
 
1504
         i2 = ifit(2,i)
 
1505
         weight = wfit(i)
 
1506
         xxyx = xxyx + weight*rion1(1,i1)*rion2(1,i2)
 
1507
         xxyy = xxyy + weight*rion1(2,i1)*rion2(1,i2)
 
1508
         xxyz = xxyz + weight*rion1(3,i1)*rion2(1,i2)
 
1509
         xyyx = xyyx + weight*rion1(1,i1)*rion2(2,i2)
 
1510
         xyyy = xyyy + weight*rion1(2,i1)*rion2(2,i2)
 
1511
         xyyz = xyyz + weight*rion1(3,i1)*rion2(2,i2)
 
1512
         xzyx = xzyx + weight*rion1(1,i1)*rion2(3,i2)
 
1513
         xzyy = xzyy + weight*rion1(2,i1)*rion2(3,i2)
 
1514
         xzyz = xzyz + weight*rion1(3,i1)*rion2(3,i2)
 
1515
      end do
 
1516
      c(1,1) = xxyx + xyyy + xzyz
 
1517
      c(1,2) = xzyy - xyyz
 
1518
      c(2,2) = xxyx - xyyy - xzyz
 
1519
      c(1,3) = xxyz - xzyx
 
1520
      c(2,3) = xxyy + xyyx
 
1521
      c(3,3) = xyyy - xzyz - xxyx
 
1522
      c(1,4) = xyyx - xxyy
 
1523
      c(2,4) = xzyx + xxyz
 
1524
      c(3,4) = xyyz + xzyy
 
1525
      c(4,4) = xzyz - xxyx - xyyy
 
1526
c
 
1527
c     diagonalize the quadratic form matrix
 
1528
c
 
1529
      call neb_jacobi4(4,4,c,d,v,temp1,temp2)
 
1530
c
 
1531
c     extract the desired quaternion
 
1532
c
 
1533
      q(1) = v(1,4)
 
1534
      q(2) = v(2,4)
 
1535
      q(3) = v(3,4)
 
1536
      q(4) = v(4,4)
 
1537
c
 
1538
c     assemble rotation matrix that superimposes the molecules
 
1539
c
 
1540
      rot(1,1) = q(1)**2 + q(2)**2 - q(3)**2 - q(4)**2
 
1541
      rot(2,1) = 2.0d0 * (q(2) * q(3) - q(1) * q(4))
 
1542
      rot(3,1) = 2.0d0 * (q(2) * q(4) + q(1) * q(3))
 
1543
      rot(1,2) = 2.0d0 * (q(3) * q(2) + q(1) * q(4))
 
1544
      rot(2,2) = q(1)**2 - q(2)**2 + q(3)**2 - q(4)**2
 
1545
      rot(3,2) = 2.0d0 * (q(3) * q(4) - q(1) * q(2))
 
1546
      rot(1,3) = 2.0d0 * (q(4) * q(2) - q(1) * q(3))
 
1547
      rot(2,3) = 2.0d0 * (q(4) * q(3) + q(1) * q(2))
 
1548
      rot(3,3) = q(1)**2 - q(2)**2 - q(3)**2 + q(4)**2
 
1549
c
 
1550
c     rotate second molecule to best fit with first molecule
 
1551
c
 
1552
      do i=1,nion
 
1553
        xrot=rion2(1,i)*rot(1,1)+rion2(2,i)*rot(1,2)+rion2(3,i)*rot(1,3)
 
1554
        drot=rion2(1,i)*rot(2,1)+rion2(2,i)*rot(2,2)+rion2(3,i)*rot(2,3)
 
1555
        zrot=rion2(1,i)*rot(3,1)+rion2(2,i)*rot(3,2)+rion2(3,i)*rot(3,3)
 
1556
        rion2(1,i) = xrot
 
1557
        rion2(2,i) = drot
 
1558
        rion2(3,i) = zrot
 
1559
      end do
 
1560
 
 
1561
      return
 
1562
      end
 
1563
 
 
1564
 
 
1565
*     ***********************************************
 
1566
*     *                                             *
 
1567
*     *             neb_jacobi4                     *
 
1568
*     *                                             *
 
1569
*     ***********************************************
 
1570
c     subroutine neb_jacobi4  --  jacobi matrix diagonalization
 
1571
c
 
1572
c
 
1573
c     This routine performs a matrix diagonalization of a real
 
1574
c     symmetric matrix by the method of Jacobi rotations
 
1575
c
 
1576
c     n    logical dimension of the matrix to be diagonalized
 
1577
c     np   physical dimension of the matrix storage area
 
1578
c     a    input with the matrix to be diagonalized; only
 
1579
c             the upper triangle and diagonal are required
 
1580
c     d    returned with the eigenvalues in ascending order
 
1581
c     v    returned with the eigenvectors of the matrix
 
1582
c     b    temporary work vector
 
1583
c     z    temporary work vector
 
1584
c
 
1585
      subroutine neb_jacobi4(n,np,a,d,v,b,z)
 
1586
      implicit none
 
1587
      integer i,j,k,ip,iq,n,np,nrot,maxrot
 
1588
      real*8  sm,tresh,s,c,t,theta,tau,h,g,p
 
1589
      real*8  a(np,np),d(np),v(np,np),b(np),z(np)
 
1590
c
 
1591
c
 
1592
c     setup and initialization
 
1593
c
 
1594
      maxrot = 100
 
1595
      nrot = 0
 
1596
      do ip = 1, n
 
1597
         do iq = 1, n
 
1598
            v(ip,iq) = 0.0d0
 
1599
         end do
 
1600
         v(ip,ip) = 1.0d0
 
1601
      end do
 
1602
      do ip = 1, n
 
1603
         b(ip) = a(ip,ip)
 
1604
         d(ip) = b(ip)
 
1605
         z(ip) = 0.0d0
 
1606
      end do
 
1607
c
 
1608
c     perform the jacobi rotations
 
1609
c
 
1610
      do i = 1, maxrot
 
1611
         sm = 0.0d0
 
1612
         do ip = 1, n-1
 
1613
            do iq = ip+1, n
 
1614
               sm = sm + abs(a(ip,iq))
 
1615
            end do
 
1616
         end do
 
1617
         if (sm .eq. 0.0d0)  goto 10
 
1618
         if (i .lt. 4) then
 
1619
            tresh = 0.2d0*sm / n**2
 
1620
         else
 
1621
            tresh = 0.0d0
 
1622
         end if
 
1623
         do ip = 1, n-1
 
1624
            do iq = ip+1, n
 
1625
               g = 100.0d0 * abs(a(ip,iq))
 
1626
               if (i.gt.4 .and. abs(d(ip))+g.eq.abs(d(ip))
 
1627
     &                    .and. abs(d(iq))+g.eq.abs(d(iq))) then
 
1628
                  a(ip,iq) = 0.0d0
 
1629
               else if (abs(a(ip,iq)) .gt. tresh) then
 
1630
                  h = d(iq) - d(ip)
 
1631
                  if (abs(h)+g .eq. abs(h)) then
 
1632
                     t = a(ip,iq) / h
 
1633
                  else
 
1634
                     theta = 0.5d0*h / a(ip,iq)
 
1635
                     t = 1.0d0 / (abs(theta)+sqrt(1.0d0+theta**2))
 
1636
                     if (theta .lt. 0.0d0)  t = -t
 
1637
                  end if
 
1638
                  c = 1.0d0 / sqrt(1.0d0+t**2)
 
1639
                  s = t * c
 
1640
                  tau = s / (1.0d0+c)
 
1641
                  h = t * a(ip,iq)
 
1642
                  z(ip) = z(ip) - h
 
1643
                  z(iq) = z(iq) + h
 
1644
                  d(ip) = d(ip) - h
 
1645
                  d(iq) = d(iq) + h
 
1646
                  a(ip,iq) = 0.0d0
 
1647
                  do j = 1, ip-1
 
1648
                     g = a(j,ip)
 
1649
                     h = a(j,iq)
 
1650
                     a(j,ip) = g - s*(h+g*tau)
 
1651
                     a(j,iq) = h + s*(g-h*tau)
 
1652
                  end do
 
1653
                  do j = ip+1, iq-1
 
1654
                     g = a(ip,j)
 
1655
                     h = a(j,iq)
 
1656
                     a(ip,j) = g - s*(h+g*tau)
 
1657
                     a(j,iq) = h + s*(g-h*tau)
 
1658
                  end do
 
1659
                  do j = iq+1, n
 
1660
                     g = a(ip,j)
 
1661
                     h = a(iq,j)
 
1662
                     a(ip,j) = g - s*(h+g*tau)
 
1663
                     a(iq,j) = h + s*(g-h*tau)
 
1664
                  end do
 
1665
                  do j = 1, n
 
1666
                     g = v(j,ip)
 
1667
                     h = v(j,iq)
 
1668
                     v(j,ip) = g - s*(h+g*tau)
 
1669
                     v(j,iq) = h + s*(g-h*tau)
 
1670
                  end do
 
1671
                  nrot = nrot + 1
 
1672
               end if
 
1673
            end do
 
1674
         end do
 
1675
         do ip = 1, n
 
1676
            b(ip) = b(ip) + z(ip)
 
1677
            d(ip) = b(ip)
 
1678
            z(ip) = 0.0d0
 
1679
         end do
 
1680
      end do
 
1681
c
 
1682
c     print warning if not converged
 
1683
c
 
1684
   10 continue
 
1685
      if (nrot.eq.maxrot) then
 
1686
         write(*,20)
 
1687
   20    format (/,' JACOBI4 -- Matrix Diagonalization not Converged')
 
1688
      end if
 
1689
c
 
1690
c     sort the eigenvalues and vectors
 
1691
c
 
1692
      do i = 1, n-1
 
1693
         k = i
 
1694
         p = d(i)
 
1695
         do j = i+1, n
 
1696
            if (d(j) .lt. p) then
 
1697
               k = j
 
1698
               p = d(j)
 
1699
            end if
 
1700
         end do
 
1701
         if (k .ne. i) then
 
1702
            d(k) = d(i)
 
1703
            d(i) = p
 
1704
            do j = 1, n
 
1705
               p = v(j,i)
 
1706
               v(j,i) = v(j,k)
 
1707
               v(j,k) = p
 
1708
            end do
 
1709
         end if
 
1710
      end do
 
1711
      return
 
1712
      end
 
1713
 
 
1714
 
 
1715
 
 
1716
 
 
1717
*     ********************************
 
1718
*     *                              *
 
1719
*     *      neb_linesearch_init     *
 
1720
*     *                              *
 
1721
*     ********************************
 
1722
      subroutine neb_linesearch_init()
 
1723
      implicit none
 
1724
 
 
1725
*     **** neblinesearch_counter common block ****
 
1726
      integer counter
 
1727
      common / neblinesearch_counter / counter
 
1728
 
 
1729
      counter = 0
 
1730
 
 
1731
      return
 
1732
      end
 
1733
 
 
1734
*     ********************************
 
1735
*     *                              *
 
1736
*     *      neb_linesearch_count    *
 
1737
*     *                              *
 
1738
*     ********************************
 
1739
      integer function neb_linesearch_count()
 
1740
      implicit none
 
1741
 
 
1742
*     **** neblinesearch_counter common block ****
 
1743
      integer counter
 
1744
      common / neblinesearch_counter / counter
 
1745
 
 
1746
      neb_linesearch_count = counter
 
1747
      return
 
1748
      end
 
1749
 
 
1750
*     ********************************
 
1751
*     *                              *
 
1752
*     *         neb_linesearch       *
 
1753
*     *                              *
 
1754
*     ********************************
 
1755
 
 
1756
      real*8 function neb_linesearch(bead_list,kbeads,t0,f0,deltat,
 
1757
     >                           tolerance,tmin,deltaE,
 
1758
     >                           stoptype)
 
1759
      implicit none
 
1760
      character*(*) bead_list
 
1761
      real*8 kbeads
 
1762
      real*8 t0,f0
 
1763
      real*8 deltat
 
1764
      real*8   tolerance
 
1765
      real*8   tmin
 
1766
      real*8   deltaE
 
1767
      integer  stoptype
 
1768
 
 
1769
 
 
1770
*     **** local variables ****
 
1771
      integer maxiter,iteration
 
1772
      parameter (maxiter=8)
 
1773
      logical secant,notfinished
 
1774
      integer indx(3)
 
1775
      real*8  t(3)
 
1776
      real*8  f(3)
 
1777
      real*8  t_last, f_last
 
1778
      real*8  t_first,f_first
 
1779
      real*8 up,down,fmin,deltaf
 
1780
 
 
1781
*     **** neblinesearch_counter common block ****
 
1782
      integer counter
 
1783
      common / neblinesearch_counter / counter
 
1784
 
 
1785
*     **** external functions ****
 
1786
      real*8   neb_line_energy
 
1787
      external neb_line_energy
 
1788
   
 
1789
      counter = counter + 1
 
1790
 
 
1791
      secant = .true.
 
1792
 
 
1793
      t(1) = t0
 
1794
      f(1)  = f0
 
1795
      t_last = t(1)
 
1796
      f_last = f(1)
 
1797
      t_first = t(1)
 
1798
      f_first = f(1)
 
1799
 
 
1800
      f(2)  =  neb_line_energy(bead_list,kbeads,t(1)+deltat,0)
 
1801
      
 
1802
      iteration = 1
 
1803
*     **** make sure that f2 < f1 ******
 
1804
      do while ((f(2).gt.f(1)).and.(iteration.le.maxiter))
 
1805
        deltat = 0.5d0*deltat
 
1806
        f(2)  =  neb_line_energy(bead_list,kbeads,t(1)+deltat,0)
 
1807
        iteration = iteration + 1
 
1808
      end do
 
1809
      t(2) = t(1) + deltat
 
1810
      t_last = t(2)
 
1811
      f_last = f(2)
 
1812
 
 
1813
 
 
1814
*     **** use secant method to generate f(3) *****
 
1815
      deltat = -f(1)*(t(2)-t(1))/(f(2)-f(1))
 
1816
      t(3)   = t(1) + deltat
 
1817
      f(3)   =  neb_line_energy(bead_list,kbeads,t(3),0)
 
1818
      iteration = iteration + 1
 
1819
      t_last = t(3)
 
1820
      f_last = f(3)
 
1821
 
 
1822
*     **** sort the function values ****
 
1823
      call neb_Order_Values(3,f,indx)
 
1824
 
 
1825
      deltaf = f(indx(2)) - f(indx(1))
 
1826
 
 
1827
 
 
1828
      if (stoptype.eq.1) then
 
1829
        notfinished = (dabs(deltaf).gt.tolerance)
 
1830
     >                .and.(iteration.le.maxiter)
 
1831
      else
 
1832
        notfinished = (dabs(f(indx(1))/f_first).gt.tolerance)
 
1833
     >                .and.(iteration.le.maxiter)
 
1834
      end if
 
1835
 
 
1836
 
 
1837
      do while (notfinished) 
 
1838
      
 
1839
 
 
1840
*       **** use secant interpolation to generate tmin ***
 
1841
        if (secant) then
 
1842
 
 
1843
          deltat = -f(indx(1))
 
1844
     >           *(t(indx(2))-t(indx(1)))
 
1845
     >           /(f(indx(2))-f(indx(1)))
 
1846
          tmin = t(indx(1)) + deltat
 
1847
          fmin  =  neb_line_energy(bead_list,kbeads,tmin,0)
 
1848
          iteration = iteration + 1
 
1849
          t_last = tmin
 
1850
          f_last = fmin
 
1851
 
 
1852
*         **** finish using secant method ****
 
1853
          if (fmin.ge.f(indx(1))) then
 
1854
            secant = .false.
 
1855
            if (fmin.lt.f(indx(3))) then
 
1856
              t(indx(3))  = tmin
 
1857
              f(indx(3))  = fmin
 
1858
              call neb_Order_Values(3,f,indx)
 
1859
            end if
 
1860
          end if
 
1861
 
 
1862
 
 
1863
        end if 
 
1864
 
 
1865
*       **** use quadradic interpolation to generate tmin ***
 
1866
        if (.not.secant) then
 
1867
          up  = (t(2)*t(2) - t(3)*t(3))*f(1)
 
1868
     >        + (t(3)*t(3) - t(1)*t(1))*f(2)
 
1869
     >        + (t(1)*t(1) - t(2)*t(2))*f(3)
 
1870
          down = (t(2) - t(3))*f(1)
 
1871
     >         + (t(3) - t(1))*f(2)
 
1872
     >         + (t(1) - t(2))*f(3)
 
1873
 
 
1874
*         **** check added by E.Apra ****
 
1875
          if(abs(down).gt.tolerance**2) then
 
1876
             tmin = 0.5d0*up/down
 
1877
             fmin  =  neb_line_energy(bead_list,kbeads,tmin,0)
 
1878
             iteration = iteration + 1
 
1879
             t_last = tmin
 
1880
             f_last = fmin
 
1881
 
 
1882
*         **** parabola fit failed - exit loop ****
 
1883
          else
 
1884
             tmin=t(indx(3))
 
1885
             fmin=f(indx(3))+tolerance
 
1886
             iteration = maxiter+1
 
1887
          endif
 
1888
 
 
1889
        end if
 
1890
 
 
1891
 
 
1892
*       **** tolerance check and replacement ****
 
1893
        if (fmin.lt.f(indx(3))) then
 
1894
           t(indx(3))  = tmin
 
1895
           f(indx(3))  = fmin
 
1896
           call neb_Order_Values(3,f,indx)
 
1897
           deltaf = f(indx(2)) - f(indx(1))
 
1898
        else
 
1899
           deltaf=0.0d0
 
1900
        end if
 
1901
 
 
1902
        if (stoptype.eq.1) then
 
1903
          notfinished = (dabs(deltaf).gt.tolerance)
 
1904
     >                .and.(iteration.le.maxiter)
 
1905
        else
 
1906
          notfinished = (dabs(f(indx(1))/f_first).gt.tolerance)
 
1907
     >                .and.(iteration.le.maxiter)
 
1908
        end if
 
1909
           
 
1910
      end do
 
1911
 
 
1912
*     **** make sure that tmin is last functions call *****
 
1913
      tmin = t(indx(1))
 
1914
      fmin = f(indx(1))
 
1915
      if (tmin.ne.t_last) 
 
1916
     >   f_last = neb_line_energy(bead_list,kbeads,tmin,0)
 
1917
 
 
1918
      deltaE = (fmin-f_first)
 
1919
  
 
1920
      neb_linesearch = fmin
 
1921
      return
 
1922
      end
 
1923
 
 
1924
*     *****************************
 
1925
*     *                           *
 
1926
*     *     neb_Order_Values      *
 
1927
*     *                           *
 
1928
*     *****************************
 
1929
*
 
1930
*   this subroutine makes f(indx(1)) < f(indx(2)) < f(indx(3)) < ....
 
1931
*   via a bubble sort
 
1932
*   
 
1933
*   Entry - n,f
 
1934
*   Exit - returns indx
 
1935
*
 
1936
 
 
1937
      subroutine neb_Order_Values(n,f,indx)
 
1938
      implicit none
 
1939
      integer n
 
1940
      real*8 f(*)
 
1941
      integer indx(*)
 
1942
 
 
1943
*     ****** local variables *****
 
1944
      integer i,j,idum
 
1945
 
 
1946
      do i=1,n
 
1947
         indx(i) = i
 
1948
      end do
 
1949
      do i=1,n-1
 
1950
        do j=i+1,n
 
1951
           if (f(indx(j)).lt.f(indx(i))) then
 
1952
              idum    = indx(i)
 
1953
              indx(i) = indx(j)
 
1954
              indx(j) = idum
 
1955
           end if
 
1956
         end do
 
1957
       end do
 
1958
            
 
1959
      return 
 
1960
      end
 
1961
 
 
1962
 
 
1963
*     *****************************
 
1964
*     *                           *
 
1965
*     *     neb_linesearch_start  *
 
1966
*     *                           *
 
1967
*     *****************************
 
1968
      subroutine neb_linesearch_start(ng_in,c0_in,s_in,
 
1969
     >                                m_in,cs_in,gs_in,
 
1970
     >                                nbeads)
 
1971
      implicit none
 
1972
      integer ng_in,m_in
 
1973
      real*8 c0_in(*)
 
1974
      real*8 s_in(*)
 
1975
      real*8 cs_in(*)
 
1976
      real*8 gs_in(*)
 
1977
      integer nbeads
 
1978
      
 
1979
#include "mafdecls.fh"
 
1980
 
 
1981
*     *** neb_linesearch ***
 
1982
      integer cs(2),gs(2),c1(2),c0(2),s(2),g1(2),t1(2),e1(2),m,ng
 
1983
      common / neblinesearch / cs,gs,c1,c0,s,g1,t1,e1,m,ng
 
1984
 
 
1985
*     **** local variables ****
 
1986
      logical value
 
1987
 
 
1988
      m = m_in + 1
 
1989
      ng = ng_in
 
1990
 
 
1991
      value  = MA_push_get(mt_dbl,m*ng,'lcs',cs(2),cs(1))
 
1992
      value  = value.and.MA_push_get(mt_dbl,m*ng,'lgs',gs(2),gs(1))
 
1993
      value  = value.and.MA_push_get(mt_dbl,ng,'lc1',c1(2),c1(1))
 
1994
      value  = value.and.MA_push_get(mt_dbl,ng,'lc0',c0(2),c0(1))
 
1995
      value  = value.and.MA_push_get(mt_dbl,ng,'ls',s(2),s(1))
 
1996
      value  = value.and.MA_push_get(mt_dbl,ng,'lg1',g1(2),g1(1))
 
1997
      value  = value.and.MA_push_get(mt_dbl,ng,'lt1',t1(2),t1(1))
 
1998
      value  = value.and.MA_push_get(mt_dbl,nbeads,'le1',e1(2),e1(1))
 
1999
      if (.not.value)
 
2000
     >   call errquit('neb_linesearch_start:increase stack memory',1,0)
 
2001
 
 
2002
      call dcopy(m_in*ng,cs_in,1,dbl_mb(cs(1)),1)
 
2003
      call dcopy(m_in*ng,gs_in,1,dbl_mb(gs(1)),1)
 
2004
      call dcopy(ng,c0_in,1,dbl_mb(c0(1)),1)
 
2005
      call dcopy(ng,s_in,1,dbl_mb(s(1)),1)
 
2006
 
 
2007
      return
 
2008
      end
 
2009
 
 
2010
*     *****************************
 
2011
*     *                           *
 
2012
*     *     neb_linesearch_end    *
 
2013
*     *                           *
 
2014
*     *****************************
 
2015
      subroutine neb_linesearch_end()
 
2016
      implicit none
 
2017
 
 
2018
#include "mafdecls.fh"
 
2019
 
 
2020
*     *** neb_linesearch ***
 
2021
      integer cs(2),gs(2),c1(2),c0(2),s(2),g1(2),t1(2),e1(2),m,ng
 
2022
      common / neblinesearch / cs,gs,c1,c0,s,g1,t1,e1,m,ng
 
2023
 
 
2024
*     **** local variables ****
 
2025
      logical value
 
2026
 
 
2027
      value =           MA_pop_stack(e1(2))
 
2028
      value = value.and.MA_pop_stack(t1(2))
 
2029
      value = value.and.MA_pop_stack(g1(2))
 
2030
      value = value.and.MA_pop_stack(s(2))
 
2031
      value = value.and.MA_pop_stack(c0(2))
 
2032
      value = value.and.MA_pop_stack(c1(2))
 
2033
      value = value.and.MA_pop_stack(gs(2))
 
2034
      value = value.and.MA_pop_stack(cs(2))
 
2035
      if (.not.value)
 
2036
     >   call errquit('neb_linesearch_end: popping stack',1,0)
 
2037
  
 
2038
      return
 
2039
      end
 
2040
 
 
2041
*     *****************************
 
2042
*     *                           *
 
2043
*     *     neb_line_energy       *
 
2044
*     *                           *
 
2045
*     *****************************
 
2046
      real*8 function neb_line_energy(bead_list,kbeads,alpha,opt)
 
2047
      implicit none
 
2048
      character*(*) bead_list
 
2049
      real*8 kbeads,alpha
 
2050
      integer opt
 
2051
 
 
2052
#include "mafdecls.fh"
 
2053
 
 
2054
*     *** neb_linesearch ***
 
2055
      integer cs(2),gs(2),c1(2),c0(2),s(2),g1(2),t1(2),e1(2),m,ng
 
2056
      common / neblinesearch / cs,gs,c1,c0,s,g1,t1,e1,m,ng
 
2057
 
 
2058
*     **** local variables ****
 
2059
      integer shift
 
2060
      real*8  ee
 
2061
 
 
2062
*     **** external functions ****
 
2063
      logical  task_gradient
 
2064
      external task_gradient
 
2065
      real*8   ddot
 
2066
      external ddot
 
2067
 
 
2068
       ee = 0.0d0
 
2069
       call neb_move(ng,
 
2070
     >               (-alpha),
 
2071
     >               dbl_mb(c0(1)),
 
2072
     >               dbl_mb(c1(1)),
 
2073
     >               dbl_mb(s(1)))
 
2074
 
 
2075
       if (opt.eq.0) then
 
2076
           call neb_coords_set(bead_list,dbl_mb(c1(1)))
 
2077
           call runmid_bead_list(bead_list,task_gradient)
 
2078
           call neb_energies_get(bead_list,dbl_mb(e1(1)))
 
2079
           call neb_gradient_get(bead_list,kbeads,
 
2080
     >                          dbl_mb(c1(1)),
 
2081
     >                          dbl_mb(e1(1)),
 
2082
     >                          dbl_mb(t1(1)),
 
2083
     >                          dbl_mb(g1(1)))
 
2084
           shift = (m-1)*ng
 
2085
           call dcopy(ng,dbl_mb(c1(1)),1,dbl_mb(cs(1)+shift),1)
 
2086
           call dcopy(ng,dbl_mb(g1(1)),1,dbl_mb(gs(1)+shift),1)
 
2087
           call neb_lmbfgs(ng,m,
 
2088
     >                     dbl_mb(cs(1)),
 
2089
     >                     dbl_mb(gs(1)),
 
2090
     >                     dbl_mb(s(1)))
 
2091
           ee = ddot(ng,dbl_mb(g1(1)),1,dbl_mb(s(1)),1)
 
2092
      end if
 
2093
 
 
2094
 
 
2095
      neb_line_energy = ee
 
2096
      return
 
2097
      end