← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwpw/pspw/qmmm/pspw_Three/pspw_Three.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/nwpw/pspw/qmmm/pspw_Three/pspw_Three.F
 
1
c
 
2
c $Id: pspw_Three.F 22780 2012-08-26 19:56:12Z bylaska $
 
3
c
 
4
 
 
5
*     *************************
 
6
*     *                       *
 
7
*     *     pspw_Three_init   *
 
8
*     *                       *
 
9
*     *************************
 
10
      subroutine pspw_Three_init(rtdb)
 
11
      implicit none
 
12
      integer rtdb
 
13
 
 
14
#include "mafdecls.fh"
 
15
#include "rtdb.fh"
 
16
#include "pspw_Three.fh"
 
17
#include "stdio.fh"
 
18
#include "errquit.fh"
 
19
 
 
20
*     **** local variables
 
21
      logical value
 
22
      integer taskid
 
23
      integer MASTER
 
24
      parameter(MASTER=0)
 
25
 
 
26
      integer i,j,l,k,nn,ntriple,shift
 
27
      integer nkatm,nkatm_qm,nkatm2
 
28
      character*50 rtdbname
 
29
 
 
30
*     **** external functions ****
 
31
      character*4 ion_atom
 
32
      character*7 c_index_name
 
33
      integer     ion_nkatm,ion_nkatm_qm
 
34
      integer     ion_katm
 
35
      external    ion_atom
 
36
      external    c_index_name
 
37
      external    ion_nkatm,ion_nkatm_qm
 
38
      external    ion_katm
 
39
      
 
40
      nkatm    = ion_nkatm()
 
41
      nkatm_qm = ion_nkatm_qm()
 
42
      nkatm2 = nkatm*nkatm
 
43
 
 
44
      value = rtdb_parallel(.true.)
 
45
 
 
46
*     **** determine the ntriples ****
 
47
      include_qm = .false.
 
48
      ntriples = 0
 
49
      do k=1,nkatm
 
50
         do j=1,nkatm
 
51
            do i=1,nkatm
 
52
 
 
53
*          **** check for standard triple potentials ****
 
54
           rtdbname='pspw_Three_ion_ion_ion_ntriple:'
 
55
     >              //ion_atom(i)//ion_atom(j)//ion_atom(k)
 
56
            if (rtdb_get(rtdb,rtdbname,mt_int,1,ntriple)) then
 
57
                ntriples = ntriples + ntriple
 
58
                if ((i.le.nkatm_qm).and.
 
59
     >              (j.le.nkatm_qm).and.
 
60
     >              (k.le.nkatm_qm)) 
 
61
     >             include_qm = .true.
 
62
            end if
 
63
         end do
 
64
      end do
 
65
 
 
66
      if (ntriples.gt.0) then
 
67
 
 
68
*        **** allocate Three parameters ****
 
69
         value = MA_alloc_get(mt_int,nkatm*nkatm*nkatm,
 
70
     >                 'ntriple_all',ntriple_all(2),ntriple_all(1))
 
71
         value = value.and.
 
72
     >        MA_alloc_get(mt_int,nkatm*nkatm*nkatm,
 
73
     >                 'triple_start',triple_start(2),triple_start(1))
 
74
         if (.not. value) 
 
75
     >   call errquit('pspw_Three_init:out of heap memory',1,MA_ERR)
 
76
         call icopy(nkatm*nkatm*nkatm,0,0,int_mb(ntriple_all(1)),1)
 
77
         call icopy(nkatm*nkatm*nkatm,0,0,int_mb(triple_start(1)),1)
 
78
 
 
79
         value = value.and.
 
80
     >        MA_alloc_get(mt_int,ntriples,
 
81
     >                 'type_all',type_all(2),type_all(1))
 
82
         value = value.and.
 
83
     >        MA_alloc_get(mt_dbl,4*ntriples,
 
84
     >                 'param_all',param_all(2),param_all(1))
 
85
         if (.not. value) 
 
86
     >   call errquit('pspw_Three_init:out of heap memory',2,MA_ERR)
 
87
 
 
88
*        **** Generate Three potential parameters ****
 
89
         nn = 0
 
90
         do k=1,nkatm
 
91
          do j=1,nkatm
 
92
            do i=1,nkatm
 
93
 
 
94
*              **** read and add standard pair potentials ****
 
95
               rtdbname='pspw_Pair_ion_ion_ion_ntriple:'
 
96
     >                  //ion_atom(i)//ion_atom(j)//ion_atom(k)
 
97
               if (.not.rtdb_get(rtdb,rtdbname,mt_int,1,ntriple))
 
98
     >            ntriple = 0
 
99
 
 
100
               shift = (k-1)*nkatm2+(j-1)*nkatm+(i-1)
 
101
               int_mb(ntriple_all(1)+shift)  = npair
 
102
               int_mb(triple_start(1)+shift) = nn
 
103
               do l=1,ntriple
 
104
                  rtdbname = 'pspw_Three_ion_ion_ion_type:'
 
105
     >                     //c_index_name(l)
 
106
     >                     //ion_atom(i)//ion_atom(j)//ion_atom(k)
 
107
                  value = value.and.
 
108
     >                    rtdb_get(rtdb,rtdbname,mt_int,1,
 
109
     >                             int_mb(type_all(1)+nn))
 
110
 
 
111
                  rtdbname = 'pspw_Three_ion_ion_ion_param:'
 
112
     >                     //c_index_name(l)
 
113
     >                     //ion_atom(i)//ion_atom(j)//ion_atom(k)
 
114
                  value = value.and.
 
115
     >                    rtdb_get(rtdb,rtdbname,mt_dbl,4,
 
116
     >                             dbl_mb(param_all(1)+4*nn))
 
117
                  nn = nn + 1
 
118
               end do
 
119
 
 
120
            end do
 
121
          end do
 
122
         end do
 
123
 
 
124
*        **** write out Three potential data ****
 
125
         call Parallel_taskid(taskid)
 
126
         if (taskid.eq.MASTER) then
 
127
             write(luout,*)
 
128
     >       'Three-Body Ion-Ion-Ion Parameters (units=a.u.):'
 
129
             if (include_qm)
 
130
     >          write(luout,*) 
 
131
     >          '- including QM/QM Three-Body interactions'
 
132
             do i=1,nkatm
 
133
             do j=i,nkatm
 
134
             do k=j,nkatm
 
135
               shift = (k-1)*nkatm2+(j-1)*nkatm+(i-1))
 
136
               ntriple = int_mb(ntriple_all(1)+shift)
 
137
               nn      = int_mb(triple_start(1)+shift)
 
138
               if (ntriple.gt.0) then
 
139
                  write(luout,'(A4,1x,A4,1x,A4)') 
 
140
     >            ion_atom(i),ion_atom(j),ion_atom(k)
 
141
 
 
142
                  do l=1,ntriple
 
143
                     if (int_mb(type_all(1)+nn+l-1).eq.1) then
 
144
 
 
145
                write(luout,'(4x,A49,E14.6,A5,E14.6,A3,E14.6,A3,E14.6)')
 
146
     >              '- Potential=A*exp(-rij/rho)-C/rij**6-D/rij**8, A:',  
 
147
     >                    dbl_mb(param_all(1)+4*(nn+l-1)),
 
148
     >                    ' rho:',  
 
149
     >                    dbl_mb(param_all(1)+4*(nn+l-1)+1),
 
150
     >                    ' C:',  
 
151
     >                    dbl_mb(param_all(1)+4*(nn+l-1)+2),
 
152
     >                    ' D:',  
 
153
     >                    dbl_mb(param_all(1)+4*(nn+l-1)+3)
 
154
 
 
155
                     end if
 
156
                  end do
 
157
               endif
 
158
             end do
 
159
             end do
 
160
             end do
 
161
             write(luout,*)
 
162
         end if
 
163
 
 
164
      end if !*** ntriples.gt.0 ****
 
165
 
 
166
      return
 
167
      end
 
168
 
 
169
*     *************************
 
170
*     *                       *
 
171
*     *     pspw_Three_end    *
 
172
*     *                       *
 
173
*     *************************
 
174
      subroutine pspw_Three_end()
 
175
      implicit none
 
176
 
 
177
#include "mafdecls.fh"
 
178
#include "pspw_Three.fh"
 
179
#include "errquit.fh"
 
180
 
 
181
      logical value
 
182
 
 
183
      if (ntriples.gt.0) then
 
184
         value =           MA_free_heap(ntriple_all(2))
 
185
         value = value.and.MA_free_heap(triple_start(2))
 
186
         value = value.and.MA_free_heap(type_all(2))
 
187
         value = value.and.MA_free_heap(param_all(2))
 
188
         if (.not.value) 
 
189
     >      call errquit('pspw_Three_end: error MA_free_heap',
 
190
     >                   0,MA_ERR)
 
191
      end if
 
192
      return
 
193
      end
 
194
 
 
195
c     *************************************
 
196
c     *                                   *
 
197
c     *       pspw_gen_Neighbor_List      *
 
198
c     *                                   *
 
199
c     *************************************
 
200
      subroutine pspw_gen_Neighbor_List(nion,rion,
 
201
     >               nshl3d,rcell,
 
202
     >               r2max,nlist_max,nlist)
 
203
      implicit none
 
204
      integer nion
 
205
      real*8 rion(3,*)
 
206
      integer nshl3d
 
207
      real*8 rcell(3,*)
 
208
      real*8 rmax
 
209
      integer nlist_max,nlist(nlist_max,*)
 
210
 
 
211
      ijkmax(1) = 0
 
212
      ijkmax(2) = 0
 
213
      ijkmax(3) = 0
 
214
      ijkmin(1) = 99999
 
215
      ijkmin(2) = 99999
 
216
      ijkmin(2) = 99999
 
217
      imax = 0
 
218
      do k=1,nion
 
219
         ijk(1,k) = nint(rion(1,k)/rmax)
 
220
         ijk(2,k) = nint(rion(2,k)/rmax)
 
221
         ijk(3,k) = nint(rion(3,k)/rmax)
 
222
 
 
223
         if (ijk(1,k).gt.ijkmax(1)) ijkmax(1) = ijk(1,k)
 
224
         if (ijk(2,k).gt.ijkmax(2)) ijkmax(2) = ijk(2,k)
 
225
         if (ijk(3,k).gt.ijkmax(3)) ijkmax(3) = ijk(3,k)
 
226
 
 
227
         if (ijk(1,k).lt.ijkmin(1)) ijkmin(1) = ijk(1,k)
 
228
         if (ijk(2,k).lt.ijkmin(2)) ijkmin(2) = ijk(2,k)
 
229
         if (ijk(3,k).lt.ijkmin(3)) ijkmin(3) = ijk(3,k)
 
230
      end do
 
231
      nijk(1) = ijkmax(1)-ijkmin(1)+1
 
232
      nijk(2) = ijkmax(2)-ijkmin(2)+1
 
233
      nijk(3) = ijkmax(3)-ijkmin(3)+1
 
234
 
 
235
      do k=ijkmin(3),ijkmax(3)
 
236
         do j=ijkmin(2),ijkmax(2)
 
237
            do i=ijkmin(1),ijkmax(1)
 
238
            end do
 
239
         end do
 
240
      end do
 
241
 
 
242
      return
 
243
      end
 
244
 
 
245
 
 
246
 
 
247
c     *************************************
 
248
c     *                                   *
 
249
c     *           pspw_Three_E            *
 
250
c     *                                   *
 
251
c     *************************************
 
252
      real*8 function pspw_Three_E(nion,nion_qm,katm,
 
253
     >                          nfrag,indx_frag_start,size_frag,kfrag,
 
254
     >                          self_interaction,lmbda,
 
255
     >                          nshl3d,rcell,
 
256
     >                          rion)
 
257
      implicit none
 
258
      integer nion,nion_qm
 
259
      integer katm(*)
 
260
      integer nfrag
 
261
      integer indx_frag_start(*),size_frag(*)
 
262
      integer kfrag(*)
 
263
      logical self_interaction(*)
 
264
      real*8  lmbda
 
265
      integer nshl3d
 
266
      real*8  rcell(nshl3d,3)
 
267
      real*8  rion(3,*)
 
268
 
 
269
#include "mafdecls.fh"
 
270
#include "pspw_Three.fh"
 
271
 
 
272
*     **** local variables ****
 
273
      integer dutask,taskid,np
 
274
      integer i,j,ii,jj,nkatm,nkatm2
 
275
      integer w1,a,k1,kk1,n1,npair,istart,k
 
276
      integer w2,b,k2,kk2
 
277
      real*8  E
 
278
 
 
279
*     **** external functions ****
 
280
      integer  ion_nkatm
 
281
      real*8   pspw_VThree_E_periodic,pspw_VThree_E_periodic_self
 
282
      real*8   pspw_VThree_E_onecell,pspw_VThree_E_periodic_image
 
283
      external ion_nkatm
 
284
      external pspw_VThree_E_periodic,pspw_VThree_E_periodic_self
 
285
      external pspw_VThree_E_onecell,pspw_VThree_E_periodic_image
 
286
      
 
287
      call nwpw_timing_start(40)
 
288
      E = 0.0d0
 
289
 
 
290
      if (ntriple.gt.0) then
 
291
 
 
292
      call Parallel_np(np)
 
293
      call Parallel_taskid(taskid)
 
294
      nkatm = ion_nkatm()
 
295
      nkatm2 = nkatm*nkatm
 
296
 
 
297
      dutask = 0
 
298
 
 
299
c     **** create neighbor list ****
 
300
      do k=1,nion
 
301
         nlist = 0
 
302
         do j=1,nion
 
303
            r2 = (rion(1,j)-rion(1,k))**2
 
304
               + (rion(2,j)-rion(2,k))**2
 
305
               + (rion(3,j)-rion(3,k))**2
 
306
            if (r2.le.r2_neighbor) then
 
307
               int_mb(list_neighbor(1)+(k-1)*nion+nlist) = j
 
308
               nlist = nlist + 1
 
309
            end if
 
310
         end do
 
311
         int_mb(nlist_neighbor(1)+k-1) = nlist
 
312
      end do
 
313
 
 
314
 
 
315
c     **** QM/QM VThree energy ****
 
316
      if (include_qm) then
 
317
         do k = 1,nion_qm-1
 
318
           if (dutask.eq.taskid) then
 
319
              kk = katm(k)
 
320
 
 
321
              nlist = int_mb(nlist_neighbor(1)+k-1)
 
322
              do r=1,nlist
 
323
                 j=int_mb(list_neighbor(1)+(k-1)*nion+r)
 
324
                 jj = katm(j)
 
325
                 do s=1,nlist
 
326
                    i=int_mb(list(1)+(k-1)*nion+s)
 
327
                    ii = katm(i)
 
328
 
 
329
                     shift = (kk-1)*nkatm2 +(jj-1)*nkatm+ii-1
 
330
                     ntriple  = int_mb(ntriple_all(1)+shift)
 
331
                     istart = int_mb(pair_start(1)+shift)
 
332
                     if (ntriple.gt.0)
 
333
     >                  E = E + pspw_VThree_E_periodic(ntriple,
 
334
     >                           int_mb(type_all(1) +istart),
 
335
     >                           dbl_mb(param_all(1)+4*istart),
 
336
     >                           rion(1,i),rion(1,j),rion(1,k)
 
337
     >                           nshl3d,rcell)
 
338
                 end do
 
339
              end do
 
340
           end if
 
341
           dutask = mod(dutask+1,np)
 
342
         end do
 
343
      end if
 
344
 
 
345
c     **** QM/MM VPower energy ****
 
346
      do j = nion_qm+1,nion
 
347
         if (dutask.eq.taskid) then
 
348
         jj = katm(j)
 
349
         do i=1,nion_qm
 
350
            ii = katm(i)
 
351
            npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
352
            istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
353
            if (npair.gt.0)
 
354
     >         E = E + pspw_VPair_E_periodic(npair,
 
355
     >                        int_mb(type_all(1) +istart),
 
356
     >                        dbl_mb(param_all(1)+4*istart),
 
357
     >                        rion(1,i),rion(1,j),
 
358
     >                        nshl3d,rcell)
 
359
         end do
 
360
         end if
 
361
         dutask = mod(dutask+1,np)
 
362
      end do
 
363
 
 
364
c     **** MM/MM VPower 1 cell energy ****
 
365
      do w1 = 1,nfrag-1
 
366
      if (dutask.eq.taskid) then
 
367
      do w2 = w1+1,nfrag
 
368
         k1 = indx_frag_start(w1)
 
369
         k2 = indx_frag_start(w2)
 
370
         kk1 = k1
 
371
         do a=1,size_frag(w1)
 
372
            kk2 = k2
 
373
            do b=1,size_frag(w2)
 
374
               ii = katm(kk1)
 
375
               jj = katm(kk2)
 
376
               npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
377
               istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
378
               if (npair.gt.0) 
 
379
     >            E = E + pspw_VPair_E_onecell(npair,
 
380
     >                        int_mb(type_all(1) +istart),
 
381
     >                        dbl_mb(param_all(1)+4*istart),
 
382
     >                        rion(1,kk1),rion(1,kk2))
 
383
               kk2 = kk2 + 1
 
384
            end do
 
385
            kk1 = kk1 + 1
 
386
         end do
 
387
      end do
 
388
      end if
 
389
      dutask = mod(dutask+1,np)
 
390
      end do
 
391
 
 
392
c     **** MM/MM VPower self energy ****
 
393
      do w1=1,nfrag
 
394
         if (self_interaction(kfrag(w1))) then
 
395
         if (dutask.eq.taskid) then
 
396
           k1 = indx_frag_start(w1)
 
397
           n1 = size_frag(w1)
 
398
           kk1 = k1
 
399
           do a=1,n1-1
 
400
             kk2 = kk1 + 1
 
401
             do b=a+1,n1
 
402
               ii = katm(kk1)
 
403
               jj = katm(kk2)
 
404
               npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
405
               istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
406
               if (npair.gt.0)
 
407
     >            E = E + pspw_VPair_E_onecell(npair,
 
408
     >                        int_mb(type_all(1) +istart),
 
409
     >                        dbl_mb(param_all(1)+4*istart),
 
410
     >                        rion(1,kk1),rion(1,kk2))
 
411
               kk2 = kk2 + 1
 
412
             end do
 
413
             kk1 = kk1 + 1
 
414
           end do
 
415
         end if
 
416
         dutask = mod(dutask+1,np)
 
417
         end if
 
418
      end do
 
419
 
 
420
c     **** MM/MM VPair self image energy ****
 
421
      if (nshl3d.gt.1) then
 
422
 
 
423
      do j = nion_qm+1,nion
 
424
       if (dutask.eq.taskid) then
 
425
         jj = katm(j)
 
426
         npair  = int_mb(npair_all(1) +(jj-1)*nkatm+jj-1)
 
427
         istart = int_mb(pair_start(1)+(jj-1)*nkatm+jj-1)
 
428
         if (npair.gt.0)
 
429
     >      E = E + pspw_VPair_E_periodic_image(npair,
 
430
     >                  int_mb(type_all(1) +istart),
 
431
     >                  dbl_mb(param_all(1)+4*istart),
 
432
     >                  nshl3d,rcell)
 
433
 
 
434
       end if
 
435
       dutask = mod(dutask+1,np)
 
436
      end do
 
437
 
 
438
 
 
439
c     **** MM/MM VPair image energy ****
 
440
      do j = (nion_qm+1),(nion-1)
 
441
         if (dutask.eq.taskid) then
 
442
         jj = katm(j)
 
443
         do i=j+1,nion
 
444
            ii = katm(i)
 
445
            npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
446
            istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
447
            if (npair.gt.0)
 
448
     >         E = E + pspw_VPair_E_periodic_self(npair,
 
449
     >                  int_mb(type_all(1) +istart),
 
450
     >                  dbl_mb(param_all(1)+4*istart),
 
451
     >                  rion(1,i),rion(1,j),
 
452
     >                  nshl3d,rcell)
 
453
         end do
 
454
         end if
 
455
         dutask = mod(dutask+1,np)
 
456
      end do
 
457
 
 
458
      end if
 
459
 
 
460
      if (np.gt.1) call Parallel_SumAll(E)
 
461
 
 
462
      end if !*** npairs.gt.0 ***
 
463
 
 
464
      call nwpw_timing_end(40)
 
465
 
 
466
      pspw_Pair_E = E
 
467
      return
 
468
      end
 
469
 
 
470
 
 
471
 
 
472
c     *********************************************
 
473
c     *                                           *
 
474
c     *              pspw_VPair_E_onecell        *
 
475
c     *                                           *
 
476
c     *********************************************
 
477
 
 
478
      real*8 function pspw_VPair_E_onecell(n,t,p,r1,r2)
 
479
      implicit none
 
480
      integer n,t(*)
 
481
      real*8 p(4,*)
 
482
      real*8 r1(3)
 
483
      real*8 r2(3)
 
484
 
 
485
*     **** local variables ****
 
486
      integer k
 
487
      real*8  dx,dy,dz,r
 
488
      real*8  E,u,u6,u12
 
489
 
 
490
      dx = r1(1) - r2(1)
 
491
      dy = r1(2) - r2(2)
 
492
      dz = r1(3) - r2(3)
 
493
      r = dsqrt(dx**2 + dy**2 + dz**2)
 
494
 
 
495
      E = 0.0d0
 
496
      do k=1,n
 
497
         if (t(k).eq.1) then
 
498
            u   = (p(2,k)/r)
 
499
            u6  = u**6
 
500
            u12 = u6**2
 
501
            E = E + 4.0d0*p(1,k)*(u12-u6)
 
502
         else if (t(k).eq.2) then
 
503
            E = E + p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6
 
504
         else if (t(k).eq.3) then
 
505
            E = E + p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6 - p(4,k)/r**8
 
506
         else if (t(k).eq.4) then
 
507
            E = E + p(1,k)*( (r-p(2,k))**p(3,k) )
 
508
         else if (t(k).eq.5) then
 
509
            E = E + p(1,k)*dexp(-r/p(2,k))
 
510
         else if (t(k).eq.6) then
 
511
            E = E + p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))**2
 
512
         endif
 
513
      end do
 
514
 
 
515
      pspw_VPair_E_onecell = E
 
516
      return
 
517
      end
 
518
 
 
519
 
 
520
c     *********************************************
 
521
c     *                                           *
 
522
c     *              pspw_VPair_E_periodic        *
 
523
c     *                                           *
 
524
c     *********************************************
 
525
 
 
526
      real*8 function pspw_VPair_E_periodic(n,t,p,r1,r2,
 
527
     >                                    nshl3d,rcell)
 
528
      implicit none
 
529
      integer n,t(*)
 
530
      real*8  p(4,*)
 
531
      real*8  r1(3)
 
532
      real*8  r2(3)
 
533
      integer nshl3d
 
534
      real*8  rcell(nshl3d,3)
 
535
 
 
536
*     **** local variables ****
 
537
      integer l,k
 
538
      real*8  dx,dy,dz
 
539
      real*8  x,y,z,r
 
540
      real*8  E,u,u6,u12
 
541
 
 
542
      E         = 0.0d0
 
543
      dx = r1(1) - r2(1)
 
544
      dy = r1(2) - r2(2)
 
545
      dz = r1(3) - r2(3)
 
546
      do l=1,nshl3d
 
547
         x = dx + rcell(l,1)
 
548
         y = dy + rcell(l,2)
 
549
         z = dz + rcell(l,3)
 
550
         r = dsqrt(x**2 + y**2 + z**2)
 
551
         do k=1,n
 
552
            if (t(k).eq.1) then
 
553
               u   = (p(2,k)/r)
 
554
               u6  = u**6
 
555
               u12 = u6**2
 
556
               E = E + 4.0d0*p(1,k)*(u12-u6)
 
557
            else if (t(k).eq.2) then
 
558
               E = E + p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6
 
559
            else if (t(k).eq.3) then
 
560
               E = E+ p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6 - p(4,k)/r**8
 
561
            else if (t(k).eq.4) then
 
562
                E = E + p(1,k)*( (r-p(2,k))**p(3,k) )
 
563
            else if (t(k).eq.5) then
 
564
                E = E + p(1,k)*dexp(-r/p(2,k))
 
565
            else if (t(k).eq.6) then
 
566
                E = E + p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))**2
 
567
            endif
 
568
         end do
 
569
      end do
 
570
      pspw_VPair_E_periodic = E
 
571
      return
 
572
      end
 
573
 
 
574
 
 
575
c     *********************************************
 
576
c     *                                           *
 
577
c     *         pspw_VPair_E_periodic_self        *
 
578
c     *                                           *
 
579
c     *********************************************
 
580
      real*8 function pspw_VPair_E_periodic_self(n,t,p,r1,r2,
 
581
     >                                         nshl3d,rcell)
 
582
      implicit none
 
583
      integer n,t(*)
 
584
      real*8 p(4,*)
 
585
      real*8 r1(3)
 
586
      real*8 r2(3)
 
587
      integer nshl3d
 
588
      real*8  rcell(nshl3d,3)
 
589
 
 
590
*     **** local variables ****
 
591
      integer l,k
 
592
      real*8  dx,dy,dz
 
593
      real*8  x,y,z,r
 
594
      real*8  E,u,u6,u12
 
595
 
 
596
      E         = 0.0d0
 
597
      dx = r1(1) - r2(1)
 
598
      dy = r1(2) - r2(2)
 
599
      dz = r1(3) - r2(3)
 
600
      do l=2,nshl3d
 
601
         x = dx + rcell(l,1)
 
602
         y = dy + rcell(l,2)
 
603
         z = dz + rcell(l,3)
 
604
         r = dsqrt(x**2 + y**2 + z**2)
 
605
         do k=1,n
 
606
            if (t(k).eq.1) then
 
607
               u   = (p(2,k)/r)
 
608
               u6  = u**6
 
609
               u12 = u6**2
 
610
               E = E + 4.0d0*p(1,k)*(u12-u6)
 
611
            else if (t(k).eq.2) then
 
612
               E = E + p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6
 
613
            else if (t(k).eq.3) then
 
614
               E = E+ p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6 - p(4,k)/r**8
 
615
            else if (t(k).eq.4) then
 
616
               E = E + p(1,k)*( (r-p(2,k))**p(3,k) )
 
617
            else if (t(k).eq.5) then
 
618
               E = E + p(1,k)*dexp(-r/p(2,k))
 
619
            else if (t(k).eq.6) then
 
620
                E = E + p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))**2
 
621
            endif
 
622
         end do
 
623
      end do
 
624
 
 
625
      pspw_VPair_E_periodic_self = E
 
626
      return
 
627
      end
 
628
 
 
629
c     *********************************************
 
630
c     *                                           *
 
631
c     *           pspw_VPair_E_periodic_image    *
 
632
c     *                                           *
 
633
c     *********************************************
 
634
      real*8 function pspw_VPair_E_periodic_image(n,t,p,nshl3d,rcell)
 
635
      implicit none
 
636
      integer n,t(*)
 
637
      real*8 p(4,*)
 
638
      integer nshl3d
 
639
      real*8  rcell(nshl3d,3)
 
640
 
 
641
*     **** local variables ****
 
642
      integer l,k
 
643
      real*8  x,y,z,r
 
644
      real*8  E,u,u6,u12
 
645
 
 
646
      E = 0.0d0
 
647
      do l=2,nshl3d
 
648
         x = rcell(l,1)
 
649
         y = rcell(l,2)
 
650
         z = rcell(l,3)
 
651
         r = dsqrt(x**2 + y**2 + z**2)
 
652
         do k=1,n
 
653
            if (t(k).eq.1) then
 
654
               u   = (p(2,k)/r)
 
655
               u6  = u**6
 
656
               u12 = u6**2
 
657
               !E = E + 4.0d0*p(1,k)*(u12-u6)
 
658
               E = E + 2.0d0*p(1,k)*(u12-u6)
 
659
            else if (t(k).eq.2) then
 
660
               !E = E + p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6
 
661
               E = E + 0.5d0*p(1,k)*dexp(-r/p(2,k)) - p(3,k)/r**6
 
662
            else if (t(k).eq.3) then
 
663
               E=E+0.5d0*p(1,k)*dexp(-r/p(2,k))-p(3,k)/r**6-p(4,k)/r**8
 
664
            else if (t(k).eq.4) then
 
665
               !E = E +p(1,k)*( (r-p(2,k))**p(3,k) )
 
666
               E = E + 0.5d0*p(1,k)*( (r-p(2,k))**p(3,k) )
 
667
            else if (t(k).eq.5) then
 
668
               !E = E + p(1,k)*dexp( r/p(2,k) )
 
669
               E = E + 0.5d0*p(1,k)*dexp(-r/p(2,k))
 
670
            else if (t(k).eq.6) then
 
671
                E = E+0.50d0*p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))**2
 
672
            endif
 
673
         end do
 
674
      end do
 
675
      pspw_VPair_E_periodic_image = E
 
676
      return
 
677
      end
 
678
 
 
679
 
 
680
 
 
681
 
 
682
c     *************************************
 
683
c     *                                   *
 
684
c     *           pspw_Pair_fion          *
 
685
c     *                                   *
 
686
c     *************************************
 
687
      subroutine pspw_Pair_fion(nion,nion_qm,katm,
 
688
     >                        nfrag,indx_frag_start,size_frag,
 
689
     >                        kfrag,
 
690
     >                        self_interaction,lmbda,
 
691
     >                        nshl3d,rcell,
 
692
     >                        rion,fion)
 
693
      implicit none
 
694
      integer nion,nion_qm
 
695
      integer katm(*)
 
696
      integer nfrag
 
697
      integer indx_frag_start(*),size_frag(*) 
 
698
      integer kfrag(*)
 
699
      logical self_interaction(*)
 
700
      real*8  lmbda
 
701
      integer nshl3d
 
702
      real*8  rcell(nshl3d,3)
 
703
      real*8  rion(3,*)
 
704
      real*8  fion(3,*)
 
705
 
 
706
#include "mafdecls.fh"
 
707
#include "pspw_Pair.fh"
 
708
 
 
709
*     **** local variables ****
 
710
      integer dutask,taskid,np
 
711
      integer i,j,ii,jj,nkatm,npair,istart
 
712
      integer w1,a,k1,kk1,n1
 
713
      integer w2,b,k2,kk2
 
714
      real*8  e1,s1,e2,s2
 
715
 
 
716
*     **** external functions ****
 
717
      integer  ion_nkatm
 
718
      external ion_nkatm
 
719
 
 
720
      call nwpw_timing_start(40)
 
721
 
 
722
      if (npairs.gt.0) then
 
723
 
 
724
      call Parallel_np(np)
 
725
      call Parallel_taskid(taskid)
 
726
      nkatm = ion_nkatm()
 
727
      dutask = 0
 
728
 
 
729
c     **** QM/QM Pair force ****
 
730
      if (include_qm) then
 
731
         do j = 1,nion_qm-1
 
732
            if (dutask.eq.taskid) then
 
733
            jj = katm(j)
 
734
            do i=j+1,nion_qm
 
735
               ii = katm(i)
 
736
               npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
737
               istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
738
               if (npair.gt.0) 
 
739
     >            call pspw_VPair_fion_periodic(npair,
 
740
     >                             int_mb(type_all(1)+istart),
 
741
     >                             dbl_mb(param_all(1)+4*istart),
 
742
     >                             rion(1,i),fion(1,i),
 
743
     >                             rion(1,j),fion(1,j),
 
744
     >                             nshl3d,rcell)
 
745
            end do
 
746
            end if
 
747
            dutask = mod(dutask+1,np)
 
748
         end do
 
749
      end if
 
750
 
 
751
c     **** QM/MM LJ energy ****
 
752
      do j = nion_qm+1,nion
 
753
         if (dutask.eq.taskid) then
 
754
         jj = katm(j)
 
755
         do i=1,nion_qm
 
756
            ii = katm(i)
 
757
            npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
758
            istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
759
            if (npair.gt.0) 
 
760
     >         call pspw_VPair_fion_periodic(npair,
 
761
     >                             int_mb(type_all(1) +istart),
 
762
     >                             dbl_mb(param_all(1)+4*istart),
 
763
     >                             rion(1,i),fion(1,i),
 
764
     >                             rion(1,j),fion(1,j),
 
765
     >                             nshl3d,rcell)
 
766
         end do
 
767
         end if
 
768
         dutask = mod(dutask+1,np)
 
769
      end do
 
770
 
 
771
c     **** MM/MM LJ 1 cell energy ****
 
772
      do w1 = 1,nfrag-1
 
773
      if (dutask.eq.taskid) then
 
774
      do w2 = w1+1,nfrag
 
775
         k1 = indx_frag_start(w1)
 
776
         k2 = indx_frag_start(w2)
 
777
         kk1 = k1
 
778
         do a=1,size_frag(w1)
 
779
            kk2 = k2
 
780
            do b=1,size_frag(w2)
 
781
               ii = katm(kk1)
 
782
               jj = katm(kk2)
 
783
               npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
784
               istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
785
               if (npair.gt.0) 
 
786
     >            call pspw_VPair_fion_onecell(npair,
 
787
     >                             int_mb(type_all(1) +istart),
 
788
     >                             dbl_mb(param_all(1)+4*istart),
 
789
     >                             rion(1,kk1),fion(1,kk1),
 
790
     >                             rion(1,kk2),fion(1,kk2))
 
791
               kk2 = kk2 + 1
 
792
            end do
 
793
            kk1 = kk1 + 1
 
794
         end do
 
795
      end do
 
796
      end if
 
797
      dutask = mod(dutask+1,np)
 
798
      end do
 
799
 
 
800
c     **** MM/MM Pair self energy ****
 
801
      do w1=1,nfrag
 
802
         if (self_interaction(kfrag(w1))) then
 
803
         if (dutask.eq.taskid) then
 
804
           k1 = indx_frag_start(w1)
 
805
           n1 = size_frag(w1)
 
806
           kk1 = k1
 
807
           do a=1,n1-1
 
808
             kk2 = kk1 + 1
 
809
             do b=a+1,n1
 
810
               ii = katm(kk1)
 
811
               jj = katm(kk2)
 
812
               npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
813
               istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
814
               if (npair.gt.0) 
 
815
     >            call pspw_VPair_fion_onecell(npair,
 
816
     >                             int_mb(type_all(1) +istart),
 
817
     >                             dbl_mb(param_all(1)+4*istart),
 
818
     >                             rion(1,kk1),fion(1,kk1),
 
819
     >                             rion(1,kk2),fion(1,kk2))
 
820
               kk2 = kk2 + 1
 
821
             end do
 
822
             kk1 = kk1 + 1
 
823
           end do
 
824
         end if
 
825
         dutask = mod(dutask+1,np)
 
826
         end if
 
827
      end do
 
828
 
 
829
      if (nshl3d.gt.1) then
 
830
 
 
831
c     **** MM/MM Pair self image energy - no force ****
 
832
c     **** MM/MM Pair image energy ****
 
833
      do j = (nion_qm+1),(nion-1)
 
834
         if (dutask.eq.taskid) then
 
835
         jj = katm(j)
 
836
         do i=j+1,nion
 
837
            ii = katm(i)
 
838
            npair  = int_mb(npair_all(1) +(jj-1)*nkatm+ii-1)
 
839
            istart = int_mb(pair_start(1)+(jj-1)*nkatm+ii-1)
 
840
            if (npair.gt.0) 
 
841
     >         call pspw_VPair_fion_periodic_self(npair,
 
842
     >                             int_mb(type_all(1) +istart),
 
843
     >                             dbl_mb(param_all(1)+4*istart),
 
844
     >                             rion(1,i),fion(1,i),
 
845
     >                             rion(1,j),fion(1,j),
 
846
     >                             nshl3d,rcell)
 
847
         end do
 
848
         end if
 
849
         dutask = mod(dutask+1,np)
 
850
      end do
 
851
 
 
852
      end if !*** nshl3d.gt.1 ***
 
853
 
 
854
      end if !*** npairs.gt.0 ***
 
855
 
 
856
 
 
857
      call nwpw_timing_end(40)
 
858
      return
 
859
      end
 
860
 
 
861
 
 
862
c     *********************************************
 
863
c     *                                           *
 
864
c     *           pspw_VPair_fion_periodic        *
 
865
c     *                                           *
 
866
c     *********************************************
 
867
      subroutine pspw_VPair_fion_periodic(n,t,p,
 
868
     >                                  r1,f1,r2,f2,
 
869
     >                                  nshl3d,rcell)
 
870
      implicit none
 
871
      integer n,t(*)
 
872
      real*8 p(4,*)
 
873
      real*8 r1(3),f1(3)
 
874
      real*8 r2(3),f2(3)
 
875
      integer nshl3d
 
876
      real*8  rcell(nshl3d,3)
 
877
 
 
878
*     **** local variables ****
 
879
      integer l,k
 
880
      real*8  dx,dy,dz
 
881
      real*8  x,y,z,r
 
882
      real*8  dVPair,u,u6,u12
 
883
 
 
884
      dx = r1(1) - r2(1)
 
885
      dy = r1(2) - r2(2)
 
886
      dz = r1(3) - r2(3)
 
887
      do l=1,nshl3d
 
888
         x = dx + rcell(l,1)
 
889
         y = dy + rcell(l,2)
 
890
         z = dz + rcell(l,3)
 
891
         r = dsqrt(x**2 + y**2 + z**2)
 
892
         do k=1,n
 
893
            if (t(k).eq.1) then
 
894
               u = (p(2,k)/r)
 
895
               u6  = u**6
 
896
               u12 = u6**2
 
897
               dVPair = -(4.0d0*p(1,k)/r)*(12.0d0*u12-6.0d0*u6)
 
898
            else if (t(k).eq.2) then
 
899
               dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
900
            else if (t(k).eq.3) then
 
901
               dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
902
     >                                               +8.0d0*p(4,k)/r**9
 
903
            else if (t(k).eq.4) then
 
904
               dVPair =  p(1,k)*p(3,k)*( (r-p(2,k))**(p(3,k)-1.0d0) )
 
905
            else if (t(k).eq.5) then
 
906
               dVPair = -p(1,k)/p(2,k)*dexp(-r/p(2,k))
 
907
            else if (t(k).eq.6) then
 
908
               dVpair = 2.0d0*p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))
 
909
     >                       *p(2,k)*dexp(-p(2,k)*(r-p(3,k)))
 
910
            else
 
911
               dVPair = 0.0d0
 
912
            endif
 
913
            f1(1) = f1(1) - (x/r)*dVPair
 
914
            f1(2) = f1(2) - (y/r)*dVPair
 
915
            f1(3) = f1(3) - (z/r)*dVPair
 
916
            f2(1) = f2(1) + (x/r)*dVPair
 
917
            f2(2) = f2(2) + (y/r)*dVPair
 
918
            f2(3) = f2(3) + (z/r)*dVPair
 
919
         end do
 
920
      end do
 
921
 
 
922
      return
 
923
      end
 
924
 
 
925
 
 
926
 
 
927
c     *********************************************
 
928
c     *                                           *
 
929
c     *      pspw_VPair_fion_periodic_self       *
 
930
c     *                                           *
 
931
c     *********************************************
 
932
      subroutine pspw_VPair_fion_periodic_self(n,t,p,
 
933
     >                                       r1,f1,r2,f2,
 
934
     >                                       nshl3d,rcell)
 
935
      implicit none
 
936
      integer n,t(*)
 
937
      real*8  p(4,*)
 
938
      real*8 r1(3),f1(3)
 
939
      real*8 r2(3),f2(3)
 
940
      integer nshl3d
 
941
      real*8  rcell(nshl3d,3)
 
942
 
 
943
*     **** local variables ****
 
944
      integer l,k
 
945
      real*8  dx,dy,dz
 
946
      real*8  x,y,z,r
 
947
      real*8  dVPair,u,u6,u12
 
948
 
 
949
      dx = r1(1) - r2(1)
 
950
      dy = r1(2) - r2(2)
 
951
      dz = r1(3) - r2(3)
 
952
      do l=2,nshl3d
 
953
         x = dx + rcell(l,1)
 
954
         y = dy + rcell(l,2)
 
955
         z = dz + rcell(l,3)
 
956
         r = dsqrt(x**2 + y**2 + z**2)
 
957
         do k=1,n
 
958
            if (t(k).eq.1) then
 
959
               u = (p(2,k)/r)
 
960
               u6  = u**6
 
961
               u12 = u6**2
 
962
               dVPair = -(4.0d0*p(1,k)/r)*(12.0d0*u12-6.0d0*u6)
 
963
            else if (t(k).eq.2) then
 
964
               dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
965
            else if (t(k).eq.3) then
 
966
               dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
967
     >                                               +8.0d0*p(4,k)/r**9
 
968
            else if (t(k).eq.4) then
 
969
               dVPair = p(1,k)*p(3,k)*( (r-p(2,k))**(p(3,k)-1.0d0) )
 
970
            else if (t(k).eq.5) then
 
971
               dVPair = -p(1,k)/p(2,k)*dexp(-r/p(2,k))
 
972
            else if (t(k).eq.6) then
 
973
               dVpair = 2.0d0*p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))
 
974
     >                       *p(2,k)*dexp(-p(2,k)*(r-p(3,k)))
 
975
            else
 
976
               dVPair = 0.0d0
 
977
            endif
 
978
            f1(1) = f1(1) - (x/r)*dVPair
 
979
            f1(2) = f1(2) - (y/r)*dVPair
 
980
            f1(3) = f1(3) - (z/r)*dVPair
 
981
            f2(1) = f2(1) + (x/r)*dVPair
 
982
            f2(2) = f2(2) + (y/r)*dVPair
 
983
            f2(3) = f2(3) + (z/r)*dVPair
 
984
         end do
 
985
      end do
 
986
 
 
987
      return
 
988
      end
 
989
 
 
990
 
 
991
c     *********************************************
 
992
c     *                                           *
 
993
c     *          pspw_VPair_fion_onecell          *
 
994
c     *                                           *
 
995
c     *********************************************
 
996
      subroutine pspw_VPair_fion_onecell(n,t,p,r1,f1,r2,f2)
 
997
      implicit none
 
998
      integer n,t(*)
 
999
      real*8 p(4,*)
 
1000
      real*8 r1(3),f1(3)
 
1001
      real*8 r2(3),f2(3)
 
1002
 
 
1003
*     **** local variables ****
 
1004
      integer k
 
1005
      real*8  x,y,z,r
 
1006
      real*8  dVPair,u,u6,u12
 
1007
 
 
1008
      x = r1(1) - r2(1)
 
1009
      y = r1(2) - r2(2)
 
1010
      z = r1(3) - r2(3)
 
1011
      r = dsqrt(x**2 + y**2 + z**2)
 
1012
      do k=1,n
 
1013
         if (t(k).eq.1) then
 
1014
            u = (p(2,k)/r)
 
1015
            u6  = u**6
 
1016
            u12 = u6**2
 
1017
            dVPair = -(4.0d0*p(1,k)/r)*(12.0d0*u12-6.0d0*u6)
 
1018
         else if (t(k).eq.2) then
 
1019
            dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
1020
         else if (t(k).eq.3) then
 
1021
            dVPair= -p(1,k)/p(2,k)*dexp(-r/p(2,k))+6.0d0*p(3,k)/r**7
 
1022
     >                                            +8.0d0*p(4,k)/r**9
 
1023
         else if (t(k).eq.4) then
 
1024
            dVPair =  p(1,k)*p(3,k)*( (r-p(2,k))**(p(3,k)-1.0d0) )
 
1025
         else if (t(k).eq.5) then
 
1026
            dVPair = -p(1,k)/p(2,k)*dexp(-r/p(2,k))
 
1027
         else if (t(k).eq.6) then
 
1028
            dVpair = 2.0d0*p(1,k)*(1.0d0-dexp(-p(2,k)*(r-p(3,k))))
 
1029
     >                    *p(2,k)*dexp(-p(2,k)*(r-p(3,k)))
 
1030
         else
 
1031
            dVpair = 0.0d0
 
1032
         endif
 
1033
         f1(1) = f1(1) - (x/r)*dVPair
 
1034
         f1(2) = f1(2) - (y/r)*dVPair
 
1035
         f1(3) = f1(3) - (z/r)*dVPair
 
1036
         f2(1) = f2(1) + (x/r)*dVPair
 
1037
         f2(2) = f2(2) + (y/r)*dVPair
 
1038
         f2(3) = f2(3) + (z/r)*dVPair
 
1039
      end do
 
1040
 
 
1041
      return
 
1042
      end
 
1043