Adam's changes
[unres.git] / source / wham / src-M-SAXS.safe / initialize_p.F
1       subroutine initialize
2
3 C Define constants and zero out tables.
4 C
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8 #ifdef MPI
9       include 'mpif.h'
10 #endif
11       include 'COMMON.IOUNITS'
12       include 'COMMON.CHAIN'
13       include 'COMMON.INTERACT'
14       include 'COMMON.GEO'
15       include 'COMMON.LOCAL'
16       include 'COMMON.TORSION'
17       include 'COMMON.FFIELD'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.MINIM' 
20       include 'COMMON.DERIV'
21       include "COMMON.WEIGHTS"
22       include "COMMON.NAMES"
23       include "COMMON.TIME1"
24       include "COMMON.TORCNSTR"
25 C
26 C The following is just to define auxiliary variables used in angle conversion
27 C
28       pi=4.0D0*datan(1.0D0)
29       dwapi=2.0D0*pi
30       dwapi3=dwapi/3.0D0
31       pipol=0.5D0*pi
32       deg2rad=pi/180.0D0
33       rad2deg=1.0D0/deg2rad
34       angmin=10.0D0*deg2rad
35 C
36 C Define I/O units.
37 C
38       inp=    1
39       iout=   2
40       ipdbin= 3
41       ipdb=   7
42       imol2=  4
43       igeom=  8
44       intin=  9
45       ithep= 11
46       irotam=12
47       itorp= 13
48       itordp= 23
49       ielep= 14
50       isidep=15 
51       isidep1=22
52       iscpp=25
53       icbase=16
54       ifourier=20
55       istat= 17
56       ientin=18
57       ientout=19
58       ibond=28
59       isccor=29
60 C
61 C WHAM files
62 C
63       ihist=30
64       iweight=31
65       izsc=32
66 C Lipidic input file for parameters range 60-79
67       iliptranpar=60
68 C
69 C Set default weights of the energy terms.
70 C
71       wlong=1.0D0
72       welec=1.0D0
73       wtor =1.0D0
74       wang =1.0D0
75       wscloc=1.0D0
76       wstrain=1.0D0
77 C
78 C Zero out tables.
79 C
80       ndih_constr=0
81       do i=1,maxres2
82         do j=1,3
83           c(j,i)=0.0D0
84           dc(j,i)=0.0D0
85         enddo
86       enddo
87       do i=1,maxres
88         do j=1,3
89           xloc(j,i)=0.0D0
90         enddo
91       enddo
92       do i=1,ntyp
93         do j=1,ntyp
94           aa_lip(i,j)=0.0D0
95           bb_lip(i,j)=0.0D0
96           aa_aq(i,j)=0.0D0
97           bb_aq(i,j)=0.0D0
98           augm(i,j)=0.0D0
99           sigma(i,j)=0.0D0
100           r0(i,j)=0.0D0
101           chi(i,j)=0.0D0
102         enddo
103         do j=1,2
104           bad(i,j)=0.0D0
105         enddo
106         chip(i)=0.0D0
107         alp(i)=0.0D0
108         sigma0(i)=0.0D0
109         sigii(i)=0.0D0
110         rr0(i)=0.0D0
111         a0thet(i)=0.0D0
112         do j=1,2
113          do ichir1=-1,1
114           do ichir2=-1,1
115           athet(j,i,ichir1,ichir2)=0.0D0
116           bthet(j,i,ichir1,ichir2)=0.0D0
117           enddo
118          enddo
119         enddo
120         do j=0,3
121           polthet(j,i)=0.0D0
122         enddo
123         do j=1,3
124           gthet(j,i)=0.0D0
125         enddo
126         theta0(i)=0.0D0
127         sig0(i)=0.0D0
128         sigc0(i)=0.0D0
129         do j=1,maxlob
130           bsc(j,i)=0.0D0
131           do k=1,3
132             censc(k,j,i)=0.0D0
133           enddo
134           do k=1,3
135             do l=1,3
136               gaussc(l,k,j,i)=0.0D0
137             enddo
138           enddo
139           nlob(i)=0
140         enddo
141       enddo
142       nlob(ntyp1)=0
143       dsc(ntyp1)=0.0D0
144       do i=-maxtor,maxtor
145         itortyp(i)=0
146        do iblock=1,2
147         do j=-maxtor,maxtor
148           do k=1,maxterm
149             v1(k,j,i,iblock)=0.0D0
150             v2(k,j,i,iblock)=0.0D0
151           enddo
152         enddo
153         enddo
154       enddo
155       do iblock=1,2
156        do i=-maxtor,maxtor
157         do j=-maxtor,maxtor
158          do k=-maxtor,maxtor
159           do l=1,maxtermd_1
160             v1c(1,l,i,j,k,iblock)=0.0D0
161             v1s(1,l,i,j,k,iblock)=0.0D0
162             v1c(2,l,i,j,k,iblock)=0.0D0
163             v1s(2,l,i,j,k,iblock)=0.0D0
164           enddo !l
165           do l=1,maxtermd_2
166            do m=1,maxtermd_2
167             v2c(m,l,i,j,k,iblock)=0.0D0
168             v2s(m,l,i,j,k,iblock)=0.0D0
169            enddo !m
170           enddo !l
171         enddo !k
172        enddo !j
173       enddo !i
174       enddo !iblock
175       do i=1,maxres
176         itype(i)=0
177         itel(i)=0
178       enddo
179 C Initialize the bridge arrays
180       ns=0
181       nss=0 
182       nhpb=0
183       do i=1,maxss
184         iss(i)=0
185       enddo
186       do i=1,maxdim
187         dhpb(i)=0.0D0
188       enddo
189       do i=1,maxres
190         ihpb(i)=0
191         jhpb(i)=0
192         dyn_ss_mask(i)=.false.
193       enddo
194 C
195 C Initialize timing.
196 C
197       call set_timers
198 C
199 C Initialize variables used in minimization.
200 C   
201 c     maxfun=5000
202 c     maxit=2000
203       maxfun=500
204       maxit=200
205       tolf=1.0D-2
206       rtolf=5.0D-4
207
208 C Initialize the variables responsible for the mode of gradient storage.
209 C
210       nfl=0
211       icg=1
212       do i=1,14
213         do j=1,14
214           if (print_order(i).eq.j) then
215             iw(print_order(i))=j
216             goto 1121
217           endif
218         enddo
219 1121    continue
220       enddo
221       calc_grad=.false.
222 C Set timers and counters for the respective routines
223       t_func = 0.0d0
224       t_grad = 0.0d0
225       t_fhel = 0.0d0
226       t_fbet = 0.0d0
227       t_ghel = 0.0d0
228       t_gbet = 0.0d0
229       t_viol = 0.0d0
230       t_gviol = 0.0d0
231       n_func = 0
232       n_grad = 0
233       n_fhel = 0
234       n_fbet = 0
235       n_ghel = 0
236       n_gbet = 0
237       n_viol = 0
238       n_gviol = 0
239       n_map = 0
240 #ifndef SPLITELE
241       nprint_ene=nprint_ene-1
242 #endif
243       return
244       end
245 c-------------------------------------------------------------------------
246       block data nazwy
247       implicit real*8 (a-h,o-z)
248       include 'DIMENSIONS'
249       include 'DIMENSIONS.ZSCOPT'
250       include 'COMMON.NAMES'
251       include 'COMMON.WEIGHTS'
252       include 'COMMON.FFIELD'
253       include 'COMMON.SHIELD'
254       data restyp /
255      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
256      & 'DSG','DGN','DSN','DTH',
257      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
258      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
259      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
260      &'AIB','ABU','D'/
261       data onelet /
262      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
263      &'a','y','w','v','l','i','f','m','c','x',
264      &'C','M','F','I','L','V','W','Y','A','G','T',
265      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
266       data potname /'LJ','LJK','BP','GB','GBV'/
267       data ename / 
268      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
269      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
270      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
271      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
272      &   "EAFM","ETHETC","ESHIELD","ESAXS"/
273       data wname /
274      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
275      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
276      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
277      &   "WLIPTRAN","WAFM","WTHETC","WSHIELD","WSAXS"/
278       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
279      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
280      &    0.0d0,0.0,0.0d0,0.0d0,0.0d0,0.0d0,0.0d0/
281       data nprint_ene /22/
282       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
283      &  16,15,17,20,21,24,22,23,26/
284       end 
285 c---------------------------------------------------------------------------
286       subroutine init_int_table
287       implicit real*8 (a-h,o-z)
288       include 'DIMENSIONS'
289       include 'DIMENSIONS.ZSCOPT'
290 #ifdef MPI
291       include 'mpif.h'
292 #endif
293 #ifdef MP
294       include 'COMMON.INFO'
295 #endif
296       include 'COMMON.CHAIN'
297       include 'COMMON.INTERACT'
298       include 'COMMON.LOCAL'
299       include 'COMMON.SBRIDGE'
300       include 'COMMON.IOUNITS'
301       include "COMMON.TORCNSTR"
302       logical scheck,lprint
303       lprint=.false.
304       do i=1,maxres
305         nint_gr(i)=0
306         nscp_gr(i)=0
307         do j=1,maxint_gr
308           istart(i,1)=0
309           iend(i,1)=0
310           ielstart(i)=0
311           ielend(i)=0
312           iscpstart(i,1)=0
313           iscpend(i,1)=0    
314         enddo
315       enddo
316       ind_scint=0
317       ind_scint_old=0
318 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
319 cd   &   (ihpb(i),jhpb(i),i=1,nss)
320       do i=nnt,nct-1
321         scheck=.false.
322         if (dyn_ss) goto 10
323         do ii=1,nss
324           if (ihpb(ii).eq.i+nres) then
325             scheck=.true.
326             jj=jhpb(ii)-nres
327             goto 10
328           endif
329         enddo
330    10   continue
331 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
332         if (scheck) then
333           if (jj.eq.i+1) then
334             nint_gr(i)=1
335             istart(i,1)=i+2
336             iend(i,1)=nct
337           else if (jj.eq.nct) then
338             nint_gr(i)=1
339             istart(i,1)=i+1
340             iend(i,1)=nct-1
341           else
342             nint_gr(i)=2
343             istart(i,1)=i+1
344             iend(i,1)=jj-1
345             istart(i,2)=jj+1
346             iend(i,2)=nct
347           endif
348         else
349           nint_gr(i)=1
350           istart(i,1)=i+1
351           iend(i,1)=nct
352           ind_scint=int_scint+nct-i
353         endif
354       enddo
355    12 continue
356       iatsc_s=nnt
357       iatsc_e=nct-1
358       if (lprint) then
359       write (iout,'(a)') 'Interaction array:'
360       do i=iatsc_s,iatsc_e
361         write (iout,'(i3,2(2x,2i3))') 
362      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
363       enddo
364       endif
365       ispp=2
366       iatel_s=nnt
367       iatel_e=nct-3
368       do i=iatel_s,iatel_e
369         ielstart(i)=i+4
370         ielend(i)=nct-1
371       enddo
372       if (lprint) then
373         write (iout,'(a)') 'Electrostatic interaction array:'
374         do i=iatel_s,iatel_e
375           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
376         enddo
377       endif ! lprint
378 c     iscp=3
379       iscp=2
380 C Partition the SC-p interaction array
381       iatscp_s=nnt
382       iatscp_e=nct-1
383       do i=nnt,nct-1
384         if (i.lt.nnt+iscp) then
385           nscp_gr(i)=1
386           iscpstart(i,1)=i+iscp
387           iscpend(i,1)=nct
388         elseif (i.gt.nct-iscp) then
389           nscp_gr(i)=1
390           iscpstart(i,1)=nnt
391           iscpend(i,1)=i-iscp
392         else
393           nscp_gr(i)=2
394           iscpstart(i,1)=nnt
395           iscpend(i,1)=i-iscp
396           iscpstart(i,2)=i+iscp
397           iscpend(i,2)=nct
398         endif 
399       enddo ! i
400       if (lprint) then
401         write (iout,'(a)') 'SC-p interaction array:'
402         do i=iatscp_s,iatscp_e
403           write (iout,'(i3,2(2x,2i3))') 
404      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
405         enddo
406       endif ! lprint
407 C Partition local interactions
408       loc_start=2
409       loc_end=nres-1
410       ithet_start=3
411       ithet_end=nres
412       iturn3_start=nnt
413       iturn3_end=nct-3
414       iturn4_start=nnt
415       iturn4_end=nct-4
416       iphi_start=nnt+3
417       iphi_end=nct
418       idihconstr_start=1
419       idihconstr_end=ndih_constr
420       ithetaconstr_start=1
421       ithetaconstr_end=ntheta_constr
422       itau_start=4
423       itau_end=nres
424       return
425       end 
426 c---------------------------------------------------------------------------
427       subroutine int_partition(int_index,lower_index,upper_index,atom,
428      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
429       implicit real*8 (a-h,o-z)
430       include 'DIMENSIONS'
431       include 'COMMON.IOUNITS'
432       integer int_index,lower_index,upper_index,atom,at_start,at_end,
433      & first_atom,last_atom,int_gr,jat_start,jat_end
434       logical lprn
435       lprn=.false.
436       if (lprn) write (iout,*) 'int_index=',int_index
437       int_index_old=int_index
438       int_index=int_index+last_atom-first_atom+1
439       if (lprn) 
440      &   write (iout,*) 'int_index=',int_index,
441      &               ' int_index_old',int_index_old,
442      &               ' lower_index=',lower_index,
443      &               ' upper_index=',upper_index,
444      &               ' atom=',atom,' first_atom=',first_atom,
445      &               ' last_atom=',last_atom
446       if (int_index.ge.lower_index) then
447         int_gr=int_gr+1
448         if (at_start.eq.0) then
449           at_start=atom
450           jat_start=first_atom-1+lower_index-int_index_old
451         else
452           jat_start=first_atom
453         endif
454         if (lprn) write (iout,*) 'jat_start',jat_start
455         if (int_index.ge.upper_index) then
456           at_end=atom
457           jat_end=first_atom-1+upper_index-int_index_old
458           return1
459         else
460           jat_end=last_atom
461         endif
462         if (lprn) write (iout,*) 'jat_end',jat_end
463       endif
464       return
465       end
466 c------------------------------------------------------------------------------
467       subroutine hpb_partition
468       implicit real*8 (a-h,o-z)
469       include 'DIMENSIONS'
470       include 'COMMON.SBRIDGE'
471       include 'COMMON.IOUNITS'
472       link_start=1
473       link_end=nhpb
474       link_start_peak=1
475       link_end_peak=npeak
476       write (iout,*) 'HPB_PARTITION',
477      &  ' nhpb',nhpb,' link_start=',link_start,
478      &  ' link_end',link_end,' link_start_peak',link_start_peak,
479      &  ' link_end_peak',link_end_peak
480       return
481       end