added source code
[unres.git] / source / unres / src_MD / src / sc_move.F
1       subroutine sc_move(n_start,n_end,n_maxtry,e_drop,
2      +     n_fun,etot)
3 c     Perform a quick search over side-chain arrangments (over
4 c     residues n_start to n_end) for a given (frozen) CA trace
5 c     Only side-chains are minimized (at most n_maxtry times each),
6 c     not CA positions
7 c     Stops if energy drops by e_drop, otherwise tries all residues
8 c     in the given range
9 c     If there is an energy drop, full minimization may be useful
10 c     n_start, n_end CAN be modified by this routine, but only if
11 c     out of bounds (n_start <= 1, n_end >= nres, n_start < n_end)
12 c     NOTE: this move should never increase the energy
13 crc      implicit none
14
15 c     Includes
16       implicit real*8 (a-h,o-z)
17       include 'DIMENSIONS'
18 #ifdef MPI
19       include 'mpif.h'
20 #endif
21       include 'COMMON.GEO'
22       include 'COMMON.VAR'
23       include 'COMMON.HEADER'
24       include 'COMMON.IOUNITS'
25       include 'COMMON.CHAIN'
26       include 'COMMON.FFIELD'
27
28 c     External functions
29       integer iran_num
30       external iran_num
31
32 c     Input arguments
33       integer n_start,n_end,n_maxtry
34       double precision e_drop
35
36 c     Output arguments
37       integer n_fun
38       double precision etot
39
40 c     Local variables
41       double precision energy(0:n_ene)
42       double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
43       double precision orig_e,cur_e
44       integer n,n_steps,n_first,n_cur,n_tot,i
45       double precision orig_w(n_ene)
46       double precision wtime
47
48
49 c     Set non side-chain weights to zero (minimization is faster)
50 c     NOTE: e(2) does not actually depend on the side-chain, only CA
51       orig_w(2)=wscp
52       orig_w(3)=welec
53       orig_w(4)=wcorr
54       orig_w(5)=wcorr5
55       orig_w(6)=wcorr6
56       orig_w(7)=wel_loc
57       orig_w(8)=wturn3
58       orig_w(9)=wturn4
59       orig_w(10)=wturn6
60       orig_w(11)=wang
61       orig_w(13)=wtor
62       orig_w(14)=wtor_d
63       orig_w(15)=wvdwpp
64
65       wscp=0.D0
66       welec=0.D0
67       wcorr=0.D0
68       wcorr5=0.D0
69       wcorr6=0.D0
70       wel_loc=0.D0
71       wturn3=0.D0
72       wturn4=0.D0
73       wturn6=0.D0
74       wang=0.D0
75       wtor=0.D0
76       wtor_d=0.D0
77       wvdwpp=0.D0
78
79 c     Make sure n_start, n_end are within proper range
80       if (n_start.lt.2) n_start=2
81       if (n_end.gt.nres-1) n_end=nres-1
82 crc      if (n_start.lt.n_end) then
83       if (n_start.gt.n_end) then
84         n_start=2
85         n_end=nres-1
86       endif
87
88 c     Save the initial values of energy and coordinates
89 cd      call chainbuild
90 cd      call etotal(energy)
91 cd      write (iout,*) 'start sc ene',energy(0)
92 cd      call enerprint(energy(0))
93 crc      etot=energy(0)
94        n_fun=0
95 crc      orig_e=etot
96 crc      cur_e=orig_e
97 crc      do i=2,nres-1
98 crc        cur_alph(i)=alph(i)
99 crc        cur_omeg(i)=omeg(i)
100 crc      enddo
101
102 ct      wtime=MPI_WTIME()
103 c     Try (one by one) all specified residues, starting from a
104 c     random position in sequence
105 c     Stop early if the energy has decreased by at least e_drop
106       n_tot=n_end-n_start+1
107       n_first=iran_num(0,n_tot-1)
108       n_steps=0
109       n=0
110 crc      do while (n.lt.n_tot .and. orig_e-etot.lt.e_drop)
111       do while (n.lt.n_tot)
112         n_cur=n_start+mod(n_first+n,n_tot)
113         call single_sc_move(n_cur,n_maxtry,e_drop,
114      +       n_steps,n_fun,etot)
115 c     If a lower energy was found, update the current structure...
116 crc        if (etot.lt.cur_e) then
117 crc          cur_e=etot
118 crc          do i=2,nres-1
119 crc            cur_alph(i)=alph(i)
120 crc            cur_omeg(i)=omeg(i)
121 crc          enddo
122 crc        else
123 c     ...else revert to the previous one
124 crc          etot=cur_e
125 crc          do i=2,nres-1
126 crc            alph(i)=cur_alph(i)
127 crc            omeg(i)=cur_omeg(i)
128 crc          enddo
129 crc        endif
130         n=n+1
131 cd
132 cd      call chainbuild
133 cd      call etotal(energy)
134 cd      print *,'running',n,energy(0)
135       enddo
136
137 cd      call chainbuild
138 cd      call etotal(energy)
139 cd      write (iout,*) 'end   sc ene',energy(0)
140
141 c     Put the original weights back to calculate the full energy
142       wscp=orig_w(2)
143       welec=orig_w(3)
144       wcorr=orig_w(4)
145       wcorr5=orig_w(5)
146       wcorr6=orig_w(6)
147       wel_loc=orig_w(7)
148       wturn3=orig_w(8)
149       wturn4=orig_w(9)
150       wturn6=orig_w(10)
151       wang=orig_w(11)
152       wtor=orig_w(13)
153       wtor_d=orig_w(14)
154       wvdwpp=orig_w(15)
155
156 crc      n_fun=n_fun+1
157 ct      write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
158       return
159       end
160
161 c-------------------------------------------------------------
162
163       subroutine single_sc_move(res_pick,n_maxtry,e_drop,
164      +     n_steps,n_fun,e_sc)
165 c     Perturb one side-chain (res_pick) and minimize the
166 c     neighbouring region, keeping all CA's and non-neighbouring
167 c     side-chains fixed
168 c     Try until e_drop energy improvement is achieved, or n_maxtry
169 c     attempts have been made
170 c     At the start, e_sc should contain the side-chain-only energy(0)
171 c     nsteps and nfun for this move are ADDED to n_steps and n_fun
172 crc      implicit none
173
174 c     Includes
175       implicit real*8 (a-h,o-z)
176       include 'DIMENSIONS'
177       include 'COMMON.VAR'
178       include 'COMMON.INTERACT'
179       include 'COMMON.CHAIN'
180       include 'COMMON.MINIM'
181       include 'COMMON.FFIELD'
182       include 'COMMON.IOUNITS'
183
184 c     External functions
185       double precision dist
186       external dist
187
188 c     Input arguments
189       integer res_pick,n_maxtry
190       double precision e_drop
191
192 c     Input/Output arguments
193       integer n_steps,n_fun
194       double precision e_sc
195
196 c     Local variables
197       logical fail
198       integer i,j
199       integer nres_moved
200       integer iretcode,loc_nfun,orig_maxfun,n_try
201       double precision sc_dist,sc_dist_cutoff
202       double precision energy(0:n_ene),orig_e,cur_e
203       double precision evdw,escloc
204       double precision cur_alph(2:nres-1),cur_omeg(2:nres-1)
205       double precision var(maxvar)
206
207       double precision orig_theta(1:nres),orig_phi(1:nres),
208      +     orig_alph(1:nres),orig_omeg(1:nres)
209
210
211 c     Define what is meant by "neighbouring side-chain"
212       sc_dist_cutoff=5.0D0
213
214 c     Don't do glycine or ends
215       i=itype(res_pick)
216       if (i.eq.10 .or. i.eq.21) return
217
218 c     Freeze everything (later will relax only selected side-chains)
219       mask_r=.true.
220       do i=1,nres
221         mask_phi(i)=0
222         mask_theta(i)=0
223         mask_side(i)=0
224       enddo
225
226 c     Find the neighbours of the side-chain to move
227 c     and save initial variables
228 crc      orig_e=e_sc
229 crc      cur_e=orig_e
230       nres_moved=0
231       do i=2,nres-1
232 c     Don't do glycine (itype(j)==10)
233         if (itype(i).ne.10) then
234           sc_dist=dist(nres+i,nres+res_pick)
235         else
236           sc_dist=sc_dist_cutoff
237         endif
238         if (sc_dist.lt.sc_dist_cutoff) then
239           nres_moved=nres_moved+1
240           mask_side(i)=1
241           cur_alph(i)=alph(i)
242           cur_omeg(i)=omeg(i)
243         endif
244       enddo
245
246       call chainbuild
247       call egb1(evdw)
248       call esc(escloc)
249       e_sc=wsc*evdw+wscloc*escloc
250 cd      call etotal(energy)
251 cd      print *,'new       ',(energy(k),k=0,n_ene)
252       orig_e=e_sc
253       cur_e=orig_e
254
255       n_try=0
256       do while (n_try.lt.n_maxtry .and. orig_e-cur_e.lt.e_drop)
257 c     Move the selected residue (don't worry if it fails)
258         call gen_side(itype(res_pick),theta(res_pick+1),
259      +       alph(res_pick),omeg(res_pick),fail)
260
261 c     Minimize the side-chains starting from the new arrangement
262         call geom_to_var(nvar,var)
263         orig_maxfun=maxfun
264         maxfun=7
265
266 crc        do i=1,nres
267 crc          orig_theta(i)=theta(i)
268 crc          orig_phi(i)=phi(i)
269 crc          orig_alph(i)=alph(i)
270 crc          orig_omeg(i)=omeg(i)
271 crc        enddo
272
273         call minimize_sc1(e_sc,var,iretcode,loc_nfun)
274         
275 cv        write(*,'(2i3,2f12.5,2i3)') 
276 cv     &       res_pick,nres_moved,orig_e,e_sc-cur_e,
277 cv     &       iretcode,loc_nfun
278
279 c$$$        if (iretcode.eq.8) then
280 c$$$          write(iout,*)'Coordinates just after code 8'
281 c$$$          call chainbuild
282 c$$$          call all_varout
283 c$$$          call flush(iout)
284 c$$$          do i=1,nres
285 c$$$            theta(i)=orig_theta(i)
286 c$$$            phi(i)=orig_phi(i)
287 c$$$            alph(i)=orig_alph(i)
288 c$$$            omeg(i)=orig_omeg(i)
289 c$$$          enddo
290 c$$$          write(iout,*)'Coordinates just before code 8'
291 c$$$          call chainbuild
292 c$$$          call all_varout
293 c$$$          call flush(iout)
294 c$$$        endif
295
296         n_fun=n_fun+loc_nfun
297         maxfun=orig_maxfun
298         call var_to_geom(nvar,var)
299
300 c     If a lower energy was found, update the current structure...
301         if (e_sc.lt.cur_e) then
302 cv              call chainbuild
303 cv              call etotal(energy)
304 cd              call egb1(evdw)
305 cd              call esc(escloc)
306 cd              e_sc1=wsc*evdw+wscloc*escloc
307 cd              print *,'     new',e_sc1,energy(0)
308 cv              print *,'new       ',energy(0)
309 cd              call enerprint(energy(0))
310           cur_e=e_sc
311           do i=2,nres-1
312             if (mask_side(i).eq.1) then
313               cur_alph(i)=alph(i)
314               cur_omeg(i)=omeg(i)
315             endif
316           enddo
317         else
318 c     ...else revert to the previous one
319           e_sc=cur_e
320           do i=2,nres-1
321             if (mask_side(i).eq.1) then
322               alph(i)=cur_alph(i)
323               omeg(i)=cur_omeg(i)
324             endif
325           enddo
326         endif
327         n_try=n_try+1
328
329       enddo
330       n_steps=n_steps+n_try
331
332 c     Reset the minimization mask_r to false
333       mask_r=.false.
334
335       return
336       end
337
338 c-------------------------------------------------------------
339
340       subroutine sc_minimize(etot,iretcode,nfun)
341 c     Minimizes side-chains only, leaving backbone frozen
342 crc      implicit none
343
344 c     Includes
345       implicit real*8 (a-h,o-z)
346       include 'DIMENSIONS'
347       include 'COMMON.VAR'
348       include 'COMMON.CHAIN'
349       include 'COMMON.FFIELD'
350
351 c     Output arguments
352       double precision etot
353       integer iretcode,nfun
354
355 c     Local variables
356       integer i
357       double precision orig_w(n_ene),energy(0:n_ene)
358       double precision var(maxvar)
359
360
361 c     Set non side-chain weights to zero (minimization is faster)
362 c     NOTE: e(2) does not actually depend on the side-chain, only CA
363       orig_w(2)=wscp
364       orig_w(3)=welec
365       orig_w(4)=wcorr
366       orig_w(5)=wcorr5
367       orig_w(6)=wcorr6
368       orig_w(7)=wel_loc
369       orig_w(8)=wturn3
370       orig_w(9)=wturn4
371       orig_w(10)=wturn6
372       orig_w(11)=wang
373       orig_w(13)=wtor
374       orig_w(14)=wtor_d
375
376       wscp=0.D0
377       welec=0.D0
378       wcorr=0.D0
379       wcorr5=0.D0
380       wcorr6=0.D0
381       wel_loc=0.D0
382       wturn3=0.D0
383       wturn4=0.D0
384       wturn6=0.D0
385       wang=0.D0
386       wtor=0.D0
387       wtor_d=0.D0
388
389 c     Prepare to freeze backbone
390       do i=1,nres
391         mask_phi(i)=0
392         mask_theta(i)=0
393         mask_side(i)=1
394       enddo
395
396 c     Minimize the side-chains
397       mask_r=.true.
398       call geom_to_var(nvar,var)
399       call minimize(etot,var,iretcode,nfun)
400       call var_to_geom(nvar,var)
401       mask_r=.false.
402
403 c     Put the original weights back and calculate the full energy
404       wscp=orig_w(2)
405       welec=orig_w(3)
406       wcorr=orig_w(4)
407       wcorr5=orig_w(5)
408       wcorr6=orig_w(6)
409       wel_loc=orig_w(7)
410       wturn3=orig_w(8)
411       wturn4=orig_w(9)
412       wturn6=orig_w(10)
413       wang=orig_w(11)
414       wtor=orig_w(13)
415       wtor_d=orig_w(14)
416
417       call chainbuild
418       call etotal(energy)
419       etot=energy(0)
420
421       return
422       end
423
424 c-------------------------------------------------------------
425       subroutine minimize_sc1(etot,x,iretcode,nfun)
426       implicit real*8 (a-h,o-z)
427       include 'DIMENSIONS'
428       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
429       include 'COMMON.IOUNITS'
430       include 'COMMON.VAR'
431       include 'COMMON.GEO'
432       include 'COMMON.MINIM'
433       common /srutu/ icall
434       dimension iv(liv)                                               
435       double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
436       double precision energia(0:n_ene)
437       external func,gradient,fdum
438       external func_restr1,grad_restr1
439       logical not_done,change,reduce 
440       common /przechowalnia/ v
441
442       call deflt(2,iv,liv,lv,v)                                         
443 * 12 means fresh start, dont call deflt                                 
444       iv(1)=12                                                          
445 * max num of fun calls                                                  
446       if (maxfun.eq.0) maxfun=500
447       iv(17)=maxfun
448 * max num of iterations                                                 
449       if (maxmin.eq.0) maxmin=1000
450       iv(18)=maxmin
451 * controls output                                                       
452       iv(19)=2                                                          
453 * selects output unit                                                   
454 c     iv(21)=iout                                                       
455       iv(21)=0
456 * 1 means to print out result                                           
457       iv(22)=0                                                          
458 * 1 means to print out summary stats                                    
459       iv(23)=0                                                          
460 * 1 means to print initial x and d                                      
461       iv(24)=0                                                          
462 * min val for v(radfac) default is 0.1                                  
463       v(24)=0.1D0                                                       
464 * max val for v(radfac) default is 4.0                                  
465       v(25)=2.0D0                                                       
466 c     v(25)=4.0D0                                                       
467 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
468 * the sumsl default is 0.1                                              
469       v(26)=0.1D0
470 * false conv if (act fnctn decrease) .lt. v(34)                         
471 * the sumsl default is 100*machep                                       
472       v(34)=v(34)/100.0D0                                               
473 * absolute convergence                                                  
474       if (tolf.eq.0.0D0) tolf=1.0D-4
475       v(31)=tolf
476 * relative convergence                                                  
477       if (rtolf.eq.0.0D0) rtolf=1.0D-4
478       v(32)=rtolf
479 * controls initial step size                                            
480        v(35)=1.0D-1                                                    
481 * large vals of d correspond to small components of step                
482       do i=1,nphi
483         d(i)=1.0D-1
484       enddo
485       do i=nphi+1,nvar
486         d(i)=1.0D-1
487       enddo
488       IF (mask_r) THEN
489        call x2xx(x,xx,nvar_restr)
490        call sumsl(nvar_restr,d,xx,func_restr1,grad_restr1,
491      &                    iv,liv,lv,v,idum,rdum,fdum)      
492        call xx2x(x,xx)
493       ELSE
494        call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)      
495       ENDIF
496       etot=v(10)                                                      
497       iretcode=iv(1)
498       nfun=iv(6)
499
500       return  
501       end  
502 ************************************************************************
503       subroutine func_restr1(n,x,nf,f,uiparm,urparm,ufparm)  
504       implicit real*8 (a-h,o-z)
505       include 'DIMENSIONS'
506       include 'COMMON.DERIV'
507       include 'COMMON.IOUNITS'
508       include 'COMMON.GEO'
509       include 'COMMON.FFIELD'
510       include 'COMMON.INTERACT'
511       include 'COMMON.TIME1'
512       common /chuju/ jjj
513       double precision energia(0:n_ene),evdw,escloc
514       integer jjj
515       double precision ufparm,e1,e2
516       external ufparm                                                   
517       integer uiparm(1)                                                 
518       real*8 urparm(1)                                                    
519       dimension x(maxvar)
520       nfl=nf
521       icg=mod(nf,2)+1
522
523 #ifdef OSF
524 c     Intercept NaNs in the coordinates, before calling etotal
525       x_sum=0.D0
526       do i=1,n
527         x_sum=x_sum+x(i)
528       enddo
529       FOUND_NAN=.false.
530       if (x_sum.ne.x_sum) then
531         write(iout,*)"   *** func_restr1 : Found NaN in coordinates"
532         f=1.0D+73
533         FOUND_NAN=.true.
534         return
535       endif
536 #endif
537
538       call var_to_geom_restr(n,x)
539       call zerograd
540       call chainbuild
541 cd    write (iout,*) 'ETOTAL called from FUNC'
542       call egb1(evdw)
543       call esc(escloc)
544       f=wsc*evdw+wscloc*escloc
545 cd      call etotal(energia(0))
546 cd      f=wsc*energia(1)+wscloc*energia(12)
547 cd      print *,f,evdw,escloc,energia(0)
548 C
549 C Sum up the components of the Cartesian gradient.
550 C
551       do i=1,nct
552         do j=1,3
553           gradx(j,i,icg)=wsc*gvdwx(j,i)
554         enddo
555       enddo
556
557       return                                                            
558       end                                                               
559 c-------------------------------------------------------
560       subroutine grad_restr1(n,x,nf,g,uiparm,urparm,ufparm)
561       implicit real*8 (a-h,o-z)
562       include 'DIMENSIONS'
563       include 'COMMON.CHAIN'
564       include 'COMMON.DERIV'
565       include 'COMMON.VAR'
566       include 'COMMON.INTERACT'
567       include 'COMMON.FFIELD'
568       include 'COMMON.IOUNITS'
569       external ufparm
570       integer uiparm(1)
571       double precision urparm(1)
572       dimension x(maxvar),g(maxvar)
573
574       icg=mod(nf,2)+1
575       if (nf-nfl+1) 20,30,40
576    20 call func_restr1(n,x,nf,f,uiparm,urparm,ufparm)
577 c     write (iout,*) 'grad 20'
578       if (nf.eq.0) return
579       goto 40
580    30 call var_to_geom_restr(n,x)
581       call chainbuild 
582 C
583 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
584 C
585    40 call cartder
586 C
587 C Convert the Cartesian gradient into internal-coordinate gradient.
588 C
589
590       ig=0
591       ind=nres-2                                                                    
592       do i=2,nres-2                
593        IF (mask_phi(i+2).eq.1) THEN                                             
594         gphii=0.0D0                                                             
595         do j=i+1,nres-1                                                         
596           ind=ind+1                                 
597           do k=1,3                                                              
598             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
599             gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
600           enddo                                                                 
601         enddo                                                                   
602         ig=ig+1
603         g(ig)=gphii
604        ELSE
605         ind=ind+nres-1-i
606        ENDIF
607       enddo                                        
608
609
610       ind=0
611       do i=1,nres-2
612        IF (mask_theta(i+2).eq.1) THEN
613         ig=ig+1
614         gthetai=0.0D0
615         do j=i+1,nres-1
616           ind=ind+1
617           do k=1,3
618             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
619             gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
620           enddo
621         enddo
622         g(ig)=gthetai
623        ELSE
624         ind=ind+nres-1-i
625        ENDIF
626       enddo
627
628       do i=2,nres-1
629         if (itype(i).ne.10) then
630          IF (mask_side(i).eq.1) THEN
631           ig=ig+1
632           galphai=0.0D0
633           do k=1,3
634             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
635           enddo
636           g(ig)=galphai
637          ENDIF
638         endif
639       enddo
640
641       
642       do i=2,nres-1
643         if (itype(i).ne.10) then
644          IF (mask_side(i).eq.1) THEN
645           ig=ig+1
646           gomegai=0.0D0
647           do k=1,3
648             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
649           enddo
650           g(ig)=gomegai
651          ENDIF
652         endif
653       enddo
654
655 C
656 C Add the components corresponding to local energy terms.
657 C
658
659       ig=0
660       igall=0
661       do i=4,nres
662         igall=igall+1
663         if (mask_phi(i).eq.1) then
664           ig=ig+1
665           g(ig)=g(ig)+gloc(igall,icg)
666         endif
667       enddo
668
669       do i=3,nres
670         igall=igall+1
671         if (mask_theta(i).eq.1) then
672           ig=ig+1
673           g(ig)=g(ig)+gloc(igall,icg)
674         endif
675       enddo
676      
677       do ij=1,2
678       do i=2,nres-1
679         if (itype(i).ne.10) then
680           igall=igall+1
681           if (mask_side(i).eq.1) then
682             ig=ig+1
683             g(ig)=g(ig)+gloc(igall,icg)
684           endif
685         endif
686       enddo
687       enddo
688
689 cd      do i=1,ig
690 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
691 cd      enddo
692       return
693       end
694 C-----------------------------------------------------------------------------
695       subroutine egb1(evdw)
696 C
697 C This subroutine calculates the interaction energy of nonbonded side chains
698 C assuming the Gay-Berne potential of interaction.
699 C
700       implicit real*8 (a-h,o-z)
701       include 'DIMENSIONS'
702       include 'COMMON.GEO'
703       include 'COMMON.VAR'
704       include 'COMMON.LOCAL'
705       include 'COMMON.CHAIN'
706       include 'COMMON.DERIV'
707       include 'COMMON.NAMES'
708       include 'COMMON.INTERACT'
709       include 'COMMON.IOUNITS'
710       include 'COMMON.CALC'
711       include 'COMMON.CONTROL'
712       logical lprn
713       evdw=0.0D0
714 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
715       evdw=0.0D0
716       lprn=.false.
717 c     if (icall.eq.0) lprn=.true.
718       ind=0
719       do i=iatsc_s,iatsc_e
720
721
722         itypi=itype(i)
723         itypi1=itype(i+1)
724         xi=c(1,nres+i)
725         yi=c(2,nres+i)
726         zi=c(3,nres+i)
727         dxi=dc_norm(1,nres+i)
728         dyi=dc_norm(2,nres+i)
729         dzi=dc_norm(3,nres+i)
730         dsci_inv=dsc_inv(itypi)
731 C
732 C Calculate SC interaction energy.
733 C
734         do iint=1,nint_gr(i)
735           do j=istart(i,iint),iend(i,iint)
736           IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
737             ind=ind+1
738             itypj=itype(j)
739             dscj_inv=dsc_inv(itypj)
740             sig0ij=sigma(itypi,itypj)
741             chi1=chi(itypi,itypj)
742             chi2=chi(itypj,itypi)
743             chi12=chi1*chi2
744             chip1=chip(itypi)
745             chip2=chip(itypj)
746             chip12=chip1*chip2
747             alf1=alp(itypi)
748             alf2=alp(itypj)
749             alf12=0.5D0*(alf1+alf2)
750 C For diagnostics only!!!
751 c           chi1=0.0D0
752 c           chi2=0.0D0
753 c           chi12=0.0D0
754 c           chip1=0.0D0
755 c           chip2=0.0D0
756 c           chip12=0.0D0
757 c           alf1=0.0D0
758 c           alf2=0.0D0
759 c           alf12=0.0D0
760             xj=c(1,nres+j)-xi
761             yj=c(2,nres+j)-yi
762             zj=c(3,nres+j)-zi
763             dxj=dc_norm(1,nres+j)
764             dyj=dc_norm(2,nres+j)
765             dzj=dc_norm(3,nres+j)
766             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
767             rij=dsqrt(rrij)
768 C Calculate angle-dependent terms of energy and contributions to their
769 C derivatives.
770             call sc_angular
771             sigsq=1.0D0/sigsq
772             sig=sig0ij*dsqrt(sigsq)
773             rij_shift=1.0D0/rij-sig+sig0ij
774 C I hate to put IF's in the loops, but here don't have another choice!!!!
775             if (rij_shift.le.0.0D0) then
776               evdw=1.0D20
777 cd              write (iout,'(2(a3,i3,2x),17(0pf7.3))')
778 cd     &        restyp(itypi),i,restyp(itypj),j,
779 cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 
780               return
781             endif
782             sigder=-sig*sigsq
783 c---------------------------------------------------------------
784             rij_shift=1.0D0/rij_shift 
785             fac=rij_shift**expon
786             e1=fac*fac*aa(itypi,itypj)
787             e2=fac*bb(itypi,itypj)
788             evdwij=eps1*eps2rt*eps3rt*(e1+e2)
789             eps2der=evdwij*eps3rt
790             eps3der=evdwij*eps2rt
791             evdwij=evdwij*eps2rt*eps3rt
792             evdw=evdw+evdwij
793             if (lprn) then
794             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
795             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
796 cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
797 cd     &        restyp(itypi),i,restyp(itypj),j,
798 cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
799 cd     &        eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
800 cd     &        om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
801 cd     &        evdwij
802             endif
803
804             if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
805      &                        'evdw',i,j,evdwij
806
807 C Calculate gradient components.
808             e1=e1*eps1*eps2rt**2*eps3rt**2
809             fac=-expon*(e1+evdwij)*rij_shift
810             sigder=fac*sigder
811             fac=rij*fac
812 C Calculate the radial part of the gradient
813             gg(1)=xj*fac
814             gg(2)=yj*fac
815             gg(3)=zj*fac
816 C Calculate angular part of the gradient.
817             call sc_grad
818           ENDIF
819           enddo      ! j
820         enddo        ! iint
821       enddo          ! i
822       end
823 C-----------------------------------------------------------------------------