2dbcc88c1293f9e71736f3f36d730e26f410a035
[unres.git] / source / unres / src_MD / src / intlocal.f
1       subroutine integral(gamma1,gamma2,gamma3,gamma4,ity1,ity2,a1,a2,
2      &  si1,si2,si3,si4,transp,q)
3       implicit none
4       integer ity1,ity2
5       integer ilam1,ilam2,ilam3,ilam4,iincr
6       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
7      &  lambda2,lambda3,lambda4
8       logical transp
9       double precision elocal,ele
10       double precision delta,delta2,sum,ene,sumene,boltz
11       double precision q,a1(2,2),a2(2,2),si1,si2,si3,si4
12       double precision conv /.01745329252d0/,pi /3.141592654d0/
13       
14       iincr=20
15       delta=iincr*conv
16       delta2=0.5d0*delta
17 cd      print *,'iincr',iincr,' delta',delta 
18 cd      write(2,*) gamma1,gamma2,ity1,ity2,a1,a2,si1,si2,si3,si4,transp
19
20 cd      do ilam1=-180,180,5
21 cd        do ilam2=-180,180,5
22 cd          lambda1=ilam1*conv+delta2
23 cd          lambda2=ilam2*conv+delta2
24 cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
25 cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
26 cd        enddo
27 cd      enddo
28 cd      stop
29
30       sum=0.0d0
31       sumene=0.0d0
32       do ilam1=-180,179,iincr
33         do ilam2=-180,179,iincr
34           do ilam3=-180,179,iincr
35             do ilam4=-180,179,iincr
36               lambda1=ilam1*conv+delta2
37               lambda2=ilam2*conv+delta2
38               lambda3=ilam3*conv+delta2
39               lambda4=ilam4*conv+delta2
40 cd              write (2,*) ilam1,ilam2,ilam3,ilam4
41 cd              write (2,*) lambda1,lambda2,lambda3,lambda4
42               ene=
43      &         -elocal(ity1,lambda1,lambda2,.false.)*
44      &          elocal(ity2,lambda3,lambda4,transp)*
45      &          ele(si1*lambda1+gamma1,si3*lambda3+gamma3,a1)*
46      &          ele(si2*lambda2+gamma2,si4*lambda4+gamma4,a2)
47 cd              write (2,*) elocal(ity1,lambda1,gamma1-pi-lambda2),
48 cd     &        elocal(ity2,lambda3,gamma2-pi-lambda4),
49 cd     &        ele(lambda1,lambda2,a1,si1,si3),
50 cd     &        ele(lambda3,lambda4,a2,si2,si4) 
51               sum=sum+ene
52             enddo
53           enddo
54         enddo
55       enddo
56       q=sum/(2*pi)**4*delta**4
57       write (2,* )'sum',sum,' q',q
58       return
59       end
60 c---------------------------------------------------------------------------
61       subroutine integral3(gamma1,gamma2,ity1,ity2,ity3,ity4,
62      &  a1,koniec,q1,q2,q3,q4)
63       implicit none
64       integer ity1,ity2,ity3,ity4
65       integer ilam1,ilam2,ilam3,ilam4,iincr
66       double precision gamma1,gamma2,gamma3,gamma4,beta,lambda1,
67      &  lambda2,lambda3,lambda4
68       logical koniec
69       double precision elocal,ele
70       double precision delta,delta2,sum1,sum2,sum3,sum4,
71      &  ene1,ene2,ene3,ene4,boltz
72       double precision q1,q2,q3,q4,a1(2,2),a2(2,2)
73       double precision conv /.01745329252d0/,pi /3.141592654d0/
74       
75       iincr=60
76       delta=iincr*conv
77       delta2=0.5d0*delta
78 cd      print *,'iincr',iincr,' delta',delta 
79       write(2,*) gamma1,gamma2,ity1,ity2,ity3,ity4,a1,koniec
80
81 cd      do ilam1=-180,180,5
82 cd        do ilam2=-180,180,5
83 cd          lambda1=ilam1*conv+delta2
84 cd          lambda2=ilam2*conv+delta2
85 cd          write(2,'(2i5,2f10.5)') ilam1,ilam2,elocal(2,lambda1,lambda2),
86 cd     &    ele(lambda1,lambda2,a1,1.0d0,1.d00)
87 cd        enddo
88 cd      enddo
89 cd      stop
90
91       sum1=0.0d0
92       sum2=0.0d0
93       sum3=0.0d0
94       sum4=0.0d0
95       do ilam1=-180,179,iincr
96         do ilam2=-180,179,iincr
97           do ilam3=-180,179,iincr
98             do ilam4=-180,179,iincr
99               lambda1=ilam1*conv+delta2
100               lambda2=ilam2*conv+delta2
101               lambda3=ilam3*conv+delta2
102               lambda4=ilam4*conv+delta2
103 cd              write (2,*) ilam1,ilam2,ilam3,ilam4
104 cd              write (2,*) lambda1,lambda2,lambda3,lambda4
105               if (.not.koniec) then
106               ene1=
107      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
108      &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
109      &          ele(lambda2,lambda4,a1)
110               else
111               ene1=
112      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
113      &          elocal(ity3,lambda3,lambda4,.false.)*
114      &          ele(lambda2,-lambda4,a1)
115               endif
116               ene2=
117      &          elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)*
118      &          elocal(ity4,lambda3,lambda4,.false.)*
119      &          ele(lambda2,lambda3,a1)
120               if (.not.koniec) then
121               ene3=
122      &          elocal(ity2,lambda1,lambda2,.false.)*
123      &          elocal(ity3,lambda3,gamma2-pi-lambda4,.false.)*
124      &          ele(lambda1,lambda4,a1)
125               else
126               ene3=
127      &          elocal(ity2,lambda1,lambda2,.false.)*
128      &          elocal(ity3,lambda3,lambda4,.false.)*
129      &          ele(lambda1,-lambda4,a1)
130               endif
131               ene4=
132      &          elocal(ity2,lambda1,lambda2,.false.)*
133      &          elocal(ity4,lambda3,lambda4,.false.)*
134      &          ele(lambda1,lambda3,a1)
135               sum1=sum1+ene1
136               sum2=sum2+ene2
137               sum3=sum3+ene3
138               sum4=sum4+ene4
139             enddo
140           enddo
141         enddo
142       enddo
143       q1=sum1/(2*pi)**4*delta**4
144       q2=sum2/(2*pi)**4*delta**4
145       q3=sum3/(2*pi)**4*delta**4
146       q4=sum4/(2*pi)**4*delta**4
147       write (2,* )'sum',sum1,sum2,sum3,sum4,' q',q1,q2,q3,q4
148       return
149       end
150 c-------------------------------------------------------------------------
151       subroutine integral5(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
152      &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4)
153       implicit none
154       integer ity1,ity2,ity3,ity4,ity5,ity6
155       integer ilam1,ilam2,ilam3,ilam4,ilam5,iincr
156       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
157      &  lambda2,lambda3,lambda4,lambda5
158       logical transp
159       double precision elocal,ele
160       double precision eloc1,eloc2,eloc3,eloc4,eloc5,eloc6,ele1,ele2
161       double precision delta,delta2,sum,ene,sumene,pom
162       double precision ene1,ene2,ene3,ene4,sum1,sum2,sum3,sum4,
163      &  a1(2,2),a2(2,2)
164       integer si1,si2,si3,si4
165       double precision conv /.01745329252d0/,pi /3.141592654d0/
166       
167       iincr=60
168       delta=iincr*conv
169       delta2=0.5d0*delta
170 cd      print *,'iincr',iincr,' delta',delta 
171 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
172 cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
173 cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
174 cd      write(2,*) 'a1=',a1
175 cd      write(2,*) 'a2=',a2
176 cd      write(2,*) si1,si2,si3,si4,transp
177
178       sum1=0.0d0
179       sum2=0.0d0
180       sum3=0.0d0
181       sum4=0.0d0
182       do ilam1=-180,179,iincr
183         do ilam2=-180,179,iincr
184           do ilam3=-180,179,iincr
185             do ilam4=-180,179,iincr
186               do ilam5=-180,179,iincr
187                 lambda1=ilam1*conv+delta2
188                 lambda2=ilam2*conv+delta2
189                 lambda3=ilam3*conv+delta2
190                 lambda4=ilam4*conv+delta2
191                 lambda5=ilam5*conv+delta2
192                 if (transp) then
193                   ele1=ele(lambda1,si4*lambda4,a1)
194                   ele2=ele(lambda2,lambda3,a2)
195                 else
196                   ele1=ele(lambda1,lambda3,a1)
197                   ele2=ele(lambda2,si4*lambda4,a2)
198                 endif
199                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
200                 eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
201                 pom=ele1*ele2*eloc2*eloc5
202                 if (si1.gt.0) then
203                   eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
204                   sum1=sum1+pom*eloc1
205                 endif
206                 eloc3=elocal(ity3,lambda2,lambda5,.false.)
207                 sum2=sum2+pom*eloc3
208                 eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
209                 sum3=sum3+pom*eloc4
210                 if (si4.gt.0) then
211                   eloc6=elocal(ity6,lambda4,lambda5,.false.)
212                   sum4=sum4+pom*eloc6
213                 endif
214               enddo
215             enddo
216           enddo
217         enddo
218       enddo
219       pom=1.0d0/(2*pi)**5*delta**5
220       ene1=sum1*pom
221       ene2=sum2*pom
222       ene3=sum3*pom
223       ene4=sum4*pom 
224 c      write (2,* )'sum',sum1,sum2,sum3,sum4,' q',ene1,ene2,ene3,ene4
225       return
226       end
227 c-------------------------------------------------------------------------
228       subroutine integral_turn6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,
229      &  ity3,ity4,ity5,ity6,a1,a2,ene_turn6)
230       implicit none
231       integer ity1,ity2,ity3,ity4,ity5,ity6
232       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
233       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
234      &  lambda2,lambda3,lambda4,lambda5,lambda6
235       logical transp
236       double precision elocal,ele
237       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
238      &  eloc61,ele1,ele2
239       double precision delta,delta2,sum,ene,sumene,pom,ene5
240       double precision ene_turn6,sum5,a1(2,2),a2(2,2)
241       double precision conv /.01745329252d0/,pi /3.141592654d0/
242       
243       iincr=60
244       delta=iincr*conv
245       delta2=0.5d0*delta
246 cd      print *,'iincr',iincr,' delta',delta 
247       write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
248      &  ' gamma3=',gamma3,' gamma4=',gamma4
249       write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
250       write(2,*) 'a1=',a1
251       write(2,*) 'a2=',a2
252
253       sum5=0.0d0
254       do ilam1=-180,179,iincr
255         do ilam2=-180,179,iincr
256           do ilam3=-180,179,iincr
257             do ilam4=-180,179,iincr
258               do ilam5=-180,179,iincr
259                 lambda1=ilam1*conv+delta2
260                 lambda2=ilam2*conv+delta2
261                 lambda3=ilam3*conv+delta2
262                 lambda4=ilam4*conv+delta2
263                 lambda5=ilam5*conv+delta2
264                 ele1=ele(lambda1,-lambda4,a1)
265                 ele2=ele(lambda2,lambda3,a2)
266                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
267                 eloc5=elocal(ity5,lambda3,lambda4,.false.)
268                 pom=ele1*ele2*eloc2*eloc5
269                 eloc3=elocal(ity3,lambda2,gamma3-pi-lambda5,.false.)
270                 eloc4=elocal(ity4,lambda5,gamma4-pi-lambda3,.false.)
271                 sum5=sum5+pom*eloc3*eloc4
272               enddo
273             enddo
274           enddo
275         enddo
276       enddo
277       pom=-1.0d0/(2*pi)**5*delta**5
278       ene_turn6=sum5*pom 
279 c      print *,'sum6',sum6,' ene6',ene6
280       return
281       end
282 c-------------------------------------------------------------------------
283       subroutine integral6(gamma1,gamma2,gamma3,gamma4,ity1,ity2,ity3,
284      &  ity4,ity5,ity6,a1,a2,si1,si2,si3,si4,transp,ene1,ene2,ene3,ene4,
285      &  ene5,ene6)
286       implicit none
287       integer ity1,ity2,ity3,ity4,ity5,ity6
288       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
289       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
290      &  lambda2,lambda3,lambda4,lambda5,lambda6
291       logical transp
292       double precision elocal,ele
293       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
294      &  eloc61,ele1,ele2
295       double precision delta,delta2,sum,ene,sumene,pom
296       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
297      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
298       integer si1,si2,si3,si4
299       double precision conv /.01745329252d0/,pi /3.141592654d0/
300       
301       iincr=60
302       delta=iincr*conv
303       delta2=0.5d0*delta
304 cd      print *,'iincr',iincr,' delta',delta 
305 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
306 cd     &  ' gamma3=',gamma3,' gamma4=',gamma4
307 cd      write(2,*) ity1,ity2,ity3,ity4,ity5,ity6
308 cd      write(2,*) 'a1=',a1
309 cd      write(2,*) 'a2=',a2
310 cd      write(2,*) si1,si2,si3,si4,transp
311
312       sum1=0.0d0
313       sum2=0.0d0
314       sum3=0.0d0
315       sum4=0.0d0
316       sum5=0.0d0
317       sum6=0.0d0
318       eloc1=0.0d0
319       eloc6=0.0d0
320       eloc61=0.0d0
321       do ilam1=-180,179,iincr
322         do ilam2=-180,179,iincr
323           do ilam3=-180,179,iincr
324             do ilam4=-180,179,iincr
325               do ilam5=-180,179,iincr
326                 do ilam6=-180,179,iincr
327                 lambda1=ilam1*conv+delta2
328                 lambda2=ilam2*conv+delta2
329                 lambda3=ilam3*conv+delta2
330                 lambda4=ilam4*conv+delta2
331                 lambda5=ilam5*conv+delta2
332                 lambda6=ilam6*conv+delta2
333                 if (transp) then
334                   ele1=ele(lambda1,si4*lambda4,a1)
335                   ele2=ele(lambda2,lambda3,a2)
336                 else
337                   ele1=ele(lambda1,lambda3,a1)
338                   ele2=ele(lambda2,si4*lambda4,a2)
339                 endif
340                 eloc2=elocal(ity2,lambda1,gamma2-pi-lambda2,.false.)
341                 eloc5=elocal(ity5,lambda3,gamma4-pi-si4*lambda4,.false.)
342                 pom=ele1*ele2*eloc2*eloc5
343                 if (si1.gt.0) then
344                   eloc1=elocal(ity1,lambda5,gamma1-pi-lambda1,.false.)
345                 endif
346                 eloc3=elocal(ity3,lambda2,lambda6,.false.)
347                 sum1=sum1+pom*eloc1*eloc3
348                 eloc4=elocal(ity4,lambda5,gamma3-pi-lambda3,.false.)
349                 if (si4.gt.0) then
350                   eloc6=elocal(ity6,lambda4,lambda6,.false.)
351                   eloc61=elocal(ity6,lambda4,lambda5,.false.)
352                 endif
353                 sum2=sum2+pom*eloc4*eloc6
354                 eloc41=elocal(ity4,lambda6,gamma3-pi-lambda3,.false.)
355                 sum3=sum3+pom*eloc1*eloc41
356                 sum4=sum4+pom*eloc1*eloc6
357                 sum5=sum5+pom*eloc3*eloc4
358                 sum6=sum6+pom*eloc3*eloc61
359                 enddo
360               enddo
361             enddo
362           enddo
363         enddo
364       enddo
365       pom=-1.0d0/(2*pi)**6*delta**6
366       ene1=sum1*pom
367       ene2=sum2*pom
368       ene3=sum3*pom
369       ene4=sum4*pom 
370       ene5=sum5*pom 
371       ene6=sum6*pom 
372 c      print *,'sum6',sum6,' ene6',ene6
373       return
374       end
375 c-------------------------------------------------------------------------
376       subroutine integral3a(gamma1,gamma2,ity1,ity2,a1,si1,ene1)
377       implicit none
378       integer ity1,ity2,ity3,ity4,ity5,ity6
379       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
380       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
381      &  lambda2,lambda3,lambda4,lambda5,lambda6
382       logical transp
383       double precision elocal,ele
384       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
385      &  eloc61,ele1,ele2
386       double precision delta,delta2,sum,ene,sumene,pom
387       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
388      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
389       integer si1,si2,si3,si4
390       double precision conv /.01745329252d0/,pi /3.141592654d0/
391       
392       iincr=60
393       delta=iincr*conv
394       delta2=0.5d0*delta
395 cd      print *,'iincr',iincr,' delta',delta 
396 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2
397 cd      write(2,*) ity1,ity2
398 cd      write(2,*) 'a1=',a1
399 cd      write(2,*) si1,
400
401       sum1=0.0d0
402       eloc1=0.0d0
403       do ilam1=-180,179,iincr
404         do ilam2=-180,179,iincr
405           do ilam3=-180,179,iincr
406             lambda1=ilam1*conv+delta2
407             lambda2=ilam2*conv+delta2
408             lambda3=ilam3*conv+delta2
409             ele1=ele(lambda1,si1*lambda3,a1)
410             eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
411             if (si1.gt.0) then
412               eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
413             else
414               eloc2=elocal(ity2,lambda2,lambda3,.false.)
415             endif
416             sum1=sum1+ele1*eloc1*eloc2
417           enddo
418         enddo
419       enddo
420       pom=1.0d0/(2*pi)**3*delta**3
421       ene1=sum1*pom
422       return
423       end
424 c-------------------------------------------------------------------------
425       subroutine integral4a(gamma1,gamma2,gamma3,ity1,ity2,ity3,a1,si1,
426      &  ene1)
427       implicit none
428       integer ity1,ity2,ity3,ity4,ity5,ity6
429       integer ilam1,ilam2,ilam3,ilam4,ilam5,ilam6,iincr
430       double precision gamma1,gamma2,gamma3,gamma4,beta,b(2,90),lambda1,
431      &  lambda2,lambda3,lambda4,lambda5,lambda6
432       logical transp
433       double precision elocal,ele
434       double precision eloc1,eloc2,eloc3,eloc4,eloc41,eloc5,eloc6,
435      &  eloc61,ele1,ele2
436       double precision delta,delta2,sum,ene,sumene,pom
437       double precision ene1,ene2,ene3,ene4,ene5,ene6,sum1,sum2,sum3,
438      &  sum4,sum5,sum6,a1(2,2),a2(2,2)
439       integer si1,si2,si3,si4
440       double precision conv /.01745329252d0/,pi /3.141592654d0/
441       
442       iincr=60
443       delta=iincr*conv
444       delta2=0.5d0*delta
445 cd      print *,'iincr',iincr,' delta',delta 
446 cd      write(2,*) 'gamma1=',gamma1,' gamma2=',gamma2,
447 cd     &  ' gamma3=',gamma3
448 cd      write(2,*) ity1,ity2,ity3
449 cd      write(2,*) 'a1=',a1
450 cd      write(2,*) 'si1=',si1
451       sum1=0.0d0
452       do ilam1=-180,179,iincr
453         do ilam2=-180,179,iincr
454           do ilam3=-180,179,iincr
455             do ilam4=-180,179,iincr
456               lambda1=ilam1*conv+delta2
457               lambda2=ilam2*conv+delta2
458               lambda3=ilam3*conv+delta2
459               lambda4=ilam4*conv+delta2
460               ele1=ele(lambda1,si1*lambda4,a1)
461               eloc1=elocal(ity1,lambda1,gamma1-pi-lambda2,.false.)
462               eloc2=elocal(ity2,lambda2,gamma2-pi-lambda3,.false.)
463               if (si1.gt.0) then
464                 eloc3=elocal(ity3,lambda3,gamma3-pi-lambda4,.false.)
465               else
466                 eloc3=elocal(ity3,lambda3,lambda4,.false.)
467               endif
468               sum1=sum1+ele1*eloc1*eloc2*eloc3
469             enddo
470           enddo
471         enddo
472       enddo
473       pom=-1.0d0/(2*pi)**4*delta**4
474       ene1=sum1*pom
475       return
476       end
477 c-------------------------------------------------------------------------
478       double precision function elocal(i,x,y,transp)
479       implicit real*8 (a-h,o-z)
480       include 'DIMENSIONS'
481       include 'COMMON.TORSION'
482       integer i
483       double precision x,y,u(2),v(2),cu(2),dv(2),ev(2) 
484       double precision scalar2
485       logical transp
486       u(1)=dcos(x)
487       u(2)=dsin(x)
488       v(1)=dcos(y)
489       v(2)=dsin(y) 
490       if (transp) then
491         call matvec2(cc(1,1,i),v,cu)
492         call matvec2(dd(1,1,i),u,dv)
493         call matvec2(ee(1,1,i),u,ev)
494         elocal=scalar2(b1(1,i),v)+scalar2(b2(1,i),u)+scalar2(cu,v)+
495      &   scalar2(dv,u)+scalar2(ev,v)
496       else 
497         call matvec2(cc(1,1,i),u,cu)
498         call matvec2(dd(1,1,i),v,dv)
499         call matvec2(ee(1,1,i),v,ev)
500         elocal=scalar2(b1(1,i),u)+scalar2(b2(1,i),v)+scalar2(cu,u)+
501      &   scalar2(dv,v)+scalar2(ev,u)
502       endif
503       return
504       end
505 c-------------------------------------------------------------------------
506       double precision function ele(x,y,a)
507       implicit none
508       double precision x,y,a(2,2),si1,si2,u(2),v(2),av(2)
509       double precision scalar2
510       u(1)=-cos(x)
511       u(2)= sin(x)
512       v(1)=-cos(y)
513       v(2)= sin(y)
514       call matvec2(a,v,av)
515       ele=scalar2(u,av) 
516       return
517       end