Merge branch 'devel' of mmka:unres into devel
[unres.git] / source / cluster / wham / src / parmread.F
1       subroutine parmread
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'sizesclu.dat'
9       include 'COMMON.IOUNITS'
10       include 'COMMON.CHAIN'
11       include 'COMMON.INTERACT'
12       include 'COMMON.GEO'
13       include 'COMMON.LOCAL'
14       include 'COMMON.TORSION'
15       include 'COMMON.FFIELD'
16       include 'COMMON.NAMES'
17       include 'COMMON.SBRIDGE'
18       include 'COMMON.SCCOR'
19       include 'COMMON.SCROT'
20       character*1 t1,t2,t3
21       character*1 onelett(4) /"G","A","P","D"/
22       logical lprint
23       dimension blower(3,3,maxlob)
24       double precision ip,mp
25 C
26 C Body
27 C
28 C Set LPRINT=.TRUE. for debugging
29       dwa16=2.0d0**(1.0d0/6.0d0)
30       lprint=.false.
31       itypro=20
32 C Assign virtual-bond length
33       vbl=3.8D0
34       vblinv=1.0D0/vbl
35       vblinv2=vblinv*vblinv
36 #ifdef CRYST_BOND
37       read (ibond,*) vbldp0,akp
38       do i=1,ntyp
39         nbondterm(i)=1
40         read (ibond,*) vbldsc0(1,i),aksc(1,i)
41         dsc(i) = vbldsc0(1,i)
42         if (i.eq.10) then
43           dsc_inv(i)=0.0D0
44         else
45           dsc_inv(i)=1.0D0/dsc(i)
46         endif
47       enddo
48 #else
49       read (ibond,*) ijunk,vbldp0,akp,rjunk
50       do i=1,ntyp
51         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
52      &   j=1,nbondterm(i))
53         dsc(i) = vbldsc0(1,i)
54         if (i.eq.10) then
55           dsc_inv(i)=0.0D0
56         else
57           dsc_inv(i)=1.0D0/dsc(i)
58         endif
59       enddo
60 #endif
61       if (lprint) then
62         write(iout,'(/a/)')"Force constants virtual bonds:"
63         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
64      &   'inertia','Pstok'
65         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
66         do i=1,ntyp
67           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
68      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
69           do j=2,nbondterm(i)
70             write (iout,'(13x,3f10.5)')
71      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
72           enddo
73         enddo
74       endif
75 #ifdef CRYST_THETA
76 C
77 C Read the parameters of the probability distribution/energy expression 
78 C of the virtual-bond valence angles theta
79 C
80       do i=1,ntyp
81         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
82         read (ithep,*) (polthet(j,i),j=0,3)
83         read (ithep,*) (gthet(j,i),j=1,3)
84         read (ithep,*) theta0(i),sig0(i),sigc0(i)
85         sigc0(i)=sigc0(i)**2
86       enddo
87       close (ithep)
88       if (lprint) then
89 c       write (iout,'(a)') 
90 c    &    'Parameters of the virtual-bond valence angles:'
91 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
92 c    & '    ATHETA0   ','         A1   ','        A2    ',
93 c    & '        B1    ','         B2   '        
94 c       do i=1,ntyp
95 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
96 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
97 c       enddo
98 c       write (iout,'(/a/9x,5a/79(1h-))') 
99 c    & 'Parameters of the expression for sigma(theta_c):',
100 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
101 c    & '     ALPH3    ','    SIGMA0C   '        
102 c       do i=1,ntyp
103 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
104 c    &      (polthet(j,i),j=0,3),sigc0(i) 
105 c       enddo
106 c       write (iout,'(/a/9x,5a/79(1h-))') 
107 c    & 'Parameters of the second gaussian:',
108 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
109 c    & '        G2    ','         G3   '        
110 c       do i=1,ntyp
111 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
112 c    &       sig0(i),(gthet(j,i),j=1,3)
113 c       enddo
114         write (iout,'(a)') 
115      &    'Parameters of the virtual-bond valence angles:'
116         write (iout,'(/a/9x,5a/79(1h-))') 
117      & 'Coefficients of expansion',
118      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
119      & '   b1*10^1    ','    b2*10^1   '        
120         do i=1,ntyp
121           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
122      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
123         enddo
124         write (iout,'(/a/9x,5a/79(1h-))') 
125      & 'Parameters of the expression for sigma(theta_c):',
126      & ' alpha0       ','  alph1       ',' alph2        ',
127      & ' alhp3        ','   sigma0c    '        
128         do i=1,ntyp
129           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
130      &      (polthet(j,i),j=0,3),sigc0(i) 
131         enddo
132         write (iout,'(/a/9x,5a/79(1h-))') 
133      & 'Parameters of the second gaussian:',
134      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
135      & '        G2    ','   G3*10^1    '        
136         do i=1,ntyp
137           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
138      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
139         enddo
140       endif
141 #else
142 C
143 C Read the parameters of Utheta determined from ab initio surfaces
144 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
145 C
146       read (ithep,*) nthetyp,ntheterm,ntheterm2,
147      &  ntheterm3,nsingle,ndouble
148       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
149       read (ithep,*) (ithetyp(i),i=1,ntyp1)
150       do i=1,maxthetyp
151         do j=1,maxthetyp
152           do k=1,maxthetyp
153             aa0thet(i,j,k)=0.0d0
154             do l=1,ntheterm
155               aathet(l,i,j,k)=0.0d0
156             enddo
157             do l=1,ntheterm2
158               do m=1,nsingle
159                 bbthet(m,l,i,j,k)=0.0d0
160                 ccthet(m,l,i,j,k)=0.0d0
161                 ddthet(m,l,i,j,k)=0.0d0
162                 eethet(m,l,i,j,k)=0.0d0
163               enddo
164             enddo
165             do l=1,ntheterm3
166               do m=1,ndouble
167                 do mm=1,ndouble
168                  ffthet(mm,m,l,i,j,k)=0.0d0
169                  ggthet(mm,m,l,i,j,k)=0.0d0
170                 enddo
171               enddo
172             enddo
173           enddo
174         enddo
175       enddo
176       do i=1,nthetyp
177         do j=1,nthetyp
178           do k=1,nthetyp
179             read (ithep,'(3a)') res1,res2,res3
180             read (ithep,*) aa0thet(i,j,k)
181             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
182             read (ithep,*)
183      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
184      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
185      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
186      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
187             read (ithep,*)
188      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
189      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
190      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
191           enddo
192         enddo
193       enddo
194 C
195 C For dummy ends assign glycine-type coefficients of theta-only terms; the
196 C coefficients of theta-and-gamma-dependent terms are zero.
197 C
198       do i=1,nthetyp
199         do j=1,nthetyp
200           do l=1,ntheterm
201             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
202             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
203           enddo
204           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
205           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
206         enddo
207         do l=1,ntheterm
208           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
209         enddo
210         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
211       enddo
212 C
213 C Control printout of the coefficients of virtual-bond-angle potentials
214 C
215       if (lprint) then
216         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
217         do i=1,nthetyp+1
218           do j=1,nthetyp+1
219             do k=1,nthetyp+1
220               write (iout,'(//4a)')
221      &         'Type ',onelett(i),onelett(j),onelett(k)
222               write (iout,'(//a,10x,a)') " l","a[l]"
223               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
224               write (iout,'(i2,1pe15.5)')
225      &           (l,aathet(l,i,j,k),l=1,ntheterm)
226             do l=1,ntheterm2
227               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
228      &          "b",l,"c",l,"d",l,"e",l
229               do m=1,nsingle
230                 write (iout,'(i2,4(1pe15.5))') m,
231      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
232      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
233               enddo
234             enddo
235             do l=1,ntheterm3
236               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
237      &          "f+",l,"f-",l,"g+",l,"g-",l
238               do m=2,ndouble
239                 do n=1,m-1
240                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
241      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
242      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
243                 enddo
244               enddo
245             enddo
246           enddo
247         enddo
248       enddo
249       call flush(iout)
250       endif
251 #endif
252
253 #ifdef CRYST_SC
254 C
255 C Read the parameters of the probability distribution/energy expression
256 C of the side chains.
257 C
258       do i=1,ntyp
259         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
260         if (i.eq.10) then
261           dsc_inv(i)=0.0D0
262         else
263           dsc_inv(i)=1.0D0/dsc(i)
264         endif
265         if (i.ne.10) then
266         do j=1,nlob(i)
267           do k=1,3
268             do l=1,3
269               blower(l,k,j)=0.0D0
270             enddo
271           enddo
272         enddo  
273         bsc(1,i)=0.0D0
274         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
275         do j=2,nlob(i)
276           read (irotam,*) bsc(j,i)
277           read (irotam,*) (censc(k,j,i),k=1,3),
278      &                                 ((blower(k,l,j),l=1,k),k=1,3)
279         enddo
280         do j=1,nlob(i)
281           do k=1,3
282             do l=1,k
283               akl=0.0D0
284               do m=1,3
285                 akl=akl+blower(k,m,j)*blower(l,m,j)
286               enddo
287               gaussc(k,l,j,i)=akl
288               gaussc(l,k,j,i)=akl
289             enddo
290           enddo 
291         enddo
292         endif
293       enddo
294       close (irotam)
295       if (lprint) then
296         write (iout,'(/a)') 'Parameters of side-chain local geometry'
297         do i=1,ntyp
298           nlobi=nlob(i)
299           if (nlobi.gt.0) then
300           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
301      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
302 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
303 c          write (iout,'(a,f10.4,4(16x,f10.4))')
304 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
305 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
306            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
307      &                             'log h',(bsc(j,i),j=1,nlobi)
308            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
309      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
310 c          write (iout,'(a)')
311 c         do j=1,nlobi
312 c           ind=0
313 c           do k=1,3
314 c             do l=1,k
315 c              ind=ind+1
316 c              blower(k,l,j)=gaussc(ind,j,i)
317 c             enddo
318 c           enddo
319 c         enddo
320           do k=1,3
321             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
322      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
323           enddo
324           endif
325         enddo
326       endif
327 #else
328 C
329 C Read scrot parameters for potentials determined from all-atom AM1 calculations
330 C added by Urszula Kozlowska 07/11/2007
331 C
332       do i=1,ntyp
333         read (irotam,*)
334        if (i.eq.10) then
335          read (irotam,*)
336        else
337          do j=1,65
338            read(irotam,*) sc_parmin(j,i)
339          enddo
340        endif
341       enddo
342 #endif
343       close(irotam)
344 #ifdef CRYST_TOR
345 C
346 C Read torsional parameters in old format
347 C
348       read (itorp,*) ntortyp,nterm_old
349       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
350       read (itorp,*) (itortyp(i),i=1,ntyp)
351       do i=1,ntortyp
352         do j=1,ntortyp
353           read (itorp,'(a)')
354           do k=1,nterm_old
355             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
356           enddo
357         enddo
358       enddo
359       close (itorp)
360       if (lprint) then
361         write (iout,'(/a/)') 'Torsional constants:'
362         do i=1,ntortyp
363           do j=1,ntortyp
364             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
365             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
366           enddo
367         enddo
368       endif
369 #else
370 C
371 C Read torsional parameters
372 C
373       read (itorp,*) ntortyp
374       read (itorp,*) (itortyp(i),i=1,ntyp)
375       write (iout,*) 'ntortyp',ntortyp
376       do i=1,ntortyp
377         do j=1,ntortyp
378           read (itorp,*) nterm(i,j),nlor(i,j)
379           v0ij=0.0d0
380           si=-1.0d0
381           do k=1,nterm(i,j)
382             read (itorp,*) kk,v1(k,i,j),v2(k,i,j) 
383             v0ij=v0ij+si*v1(k,i,j)
384             si=-si
385           enddo
386           do k=1,nlor(i,j)
387             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) 
388             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
389           enddo
390           v0(i,j)=v0ij
391         enddo
392       enddo
393       close (itorp)
394       if (lprint) then
395         write (iout,'(/a/)') 'Torsional constants:'
396         do i=1,ntortyp
397           do j=1,ntortyp
398             write (iout,*) 'ityp',i,' jtyp',j
399             write (iout,*) 'Fourier constants'
400             do k=1,nterm(i,j)
401               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
402             enddo
403             write (iout,*) 'Lorenz constants'
404             do k=1,nlor(i,j)
405               write (iout,'(3(1pe15.5))') 
406      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
407             enddo
408           enddo
409         enddo
410       endif
411 C
412 C 6/23/01 Read parameters for double torsionals
413 C
414       do i=1,ntortyp
415         do j=1,ntortyp
416           do k=1,ntortyp
417             read (itordp,'(3a1)') t1,t2,t3
418             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
419      &        .or. t3.ne.onelett(k)) then
420               write (iout,*) "Error in double torsional parameter file",
421      &         i,j,k,t1,t2,t3
422                stop "Error in double torsional parameter file"
423             endif
424             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
425             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
426             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
427             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
428             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
429             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
430      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
431           enddo
432         enddo
433       enddo
434       if (lprint) then
435       write (iout,*) 
436       write (iout,*) 'Constants for double torsionals'
437       do i=1,ntortyp
438         do j=1,ntortyp 
439           do k=1,ntortyp
440             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
441      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
442             write (iout,*)
443             write (iout,*) 'Single angles:'
444             do l=1,ntermd_1(i,j,k)
445               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
446      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
447      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
448             enddo
449             write (iout,*)
450             write (iout,*) 'Pairs of angles:'
451             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
452             do l=1,ntermd_2(i,j,k)
453               write (iout,'(i5,20f10.5)') 
454      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
455             enddo
456             write (iout,*)
457             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
458             do l=1,ntermd_2(i,j,k)
459               write (iout,'(i5,20f10.5)') 
460      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
461             enddo
462             write (iout,*)
463           enddo
464         enddo
465       enddo
466       endif
467 #endif
468 C Read of Side-chain backbone correlation parameters
469 C Modified 11 May 2012 by Adasko
470 CCC
471 C
472       read (isccor,*) nsccortyp
473       read (isccor,*) (isccortyp(i),i=1,ntyp)
474 c      write (iout,*) 'ntortyp',ntortyp
475       maxinter=3
476 cc maxinter is maximum interaction sites
477       do l=1,maxinter
478       do i=1,nsccortyp
479         do j=1,nsccortyp
480           read (isccor,*)
481      &nterm_sccor(i,j),nlor_sccor(i,j)
482           v0ijsccor=0.0d0
483           si=-1.0d0
484
485           do k=1,nterm_sccor(i,j)
486             read (isccor,*) kk,v1sccor(k,l,i,j)
487      &    ,v2sccor(k,l,i,j)
488             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
489             si=-si
490           enddo
491           do k=1,nlor_sccor(i,j)
492             read (isccor,*) kk,vlor1sccor(k,i,j),
493      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
494             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
495      &(1+vlor3sccor(k,i,j)**2)
496           enddo
497           v0sccor(i,j)=v0ijsccor
498         enddo
499       enddo
500       enddo
501       close (isccor)
502
503       if (lprint) then
504         write (iout,'(/a/)') 'Torsional constants:'
505         do i=1,nsccortyp
506           do j=1,nsccortyp
507             write (iout,*) 'ityp',i,' jtyp',j
508             write (iout,*) 'Fourier constants'
509             do k=1,nterm_sccor(i,j)
510               write (iout,'(2(1pe15.5))')
511      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
512             enddo
513             write (iout,*) 'Lorenz constants'
514             do k=1,nlor_sccor(i,j)
515               write (iout,'(3(1pe15.5))')
516      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
517             enddo
518           enddo
519         enddo
520       endif
521
522 C
523 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
524 C         interaction energy of the Gly, Ala, and Pro prototypes.
525 C
526       read (ifourier,*) nloctyp
527       do i=1,nloctyp
528         read (ifourier,*)
529         read (ifourier,*) (b(ii,i),ii=1,13)
530         if (lprint) then
531         write (iout,*) 'Type',i
532         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
533         endif
534         B1(1,i)  = b(3,i)
535         B1(2,i)  = b(5,i)
536         B1tilde(1,i) = b(3,i)
537         B1tilde(2,i) =-b(5,i) 
538         B2(1,i)  = b(2,i)
539         B2(2,i)  = b(4,i)
540         CC(1,1,i)= b(7,i)
541         CC(2,2,i)=-b(7,i)
542         CC(2,1,i)= b(9,i)
543         CC(1,2,i)= b(9,i)
544         Ctilde(1,1,i)=b(7,i)
545         Ctilde(1,2,i)=b(9,i)
546         Ctilde(2,1,i)=-b(9,i)
547         Ctilde(2,2,i)=b(7,i)
548         DD(1,1,i)= b(6,i)
549         DD(2,2,i)=-b(6,i)
550         DD(2,1,i)= b(8,i)
551         DD(1,2,i)= b(8,i)
552         Dtilde(1,1,i)=b(6,i)
553         Dtilde(1,2,i)=b(8,i)
554         Dtilde(2,1,i)=-b(8,i)
555         Dtilde(2,2,i)=b(6,i)
556         EE(1,1,i)= b(10,i)+b(11,i)
557         EE(2,2,i)=-b(10,i)+b(11,i)
558         EE(2,1,i)= b(12,i)-b(13,i)
559         EE(1,2,i)= b(12,i)+b(13,i)
560       enddo
561       if (lprint) then
562       do i=1,nloctyp
563         write (iout,*) 'Type',i
564         write (iout,*) 'B1'
565 c        write (iout,'(f10.5)') B1(:,i)
566         write(iout,*) B1(1,i),B1(2,i)
567         write (iout,*) 'B2'
568 c        write (iout,'(f10.5)') B2(:,i)
569         write(iout,*) B2(1,i),B2(2,i)
570         write (iout,*) 'CC'
571         do j=1,2
572           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
573         enddo
574         write(iout,*) 'DD'
575         do j=1,2
576           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
577         enddo
578         write(iout,*) 'EE'
579         do j=1,2
580           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
581         enddo
582       enddo
583       endif
584
585 C Read electrostatic-interaction parameters
586 C
587       if (lprint) then
588         write (iout,'(/a)') 'Electrostatic interaction constants:'
589         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
590      &            'IT','JT','APP','BPP','AEL6','AEL3'
591       endif
592       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
593       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
594       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
595       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
596       close (ielep)
597       do i=1,2
598         do j=1,2
599         rri=rpp(i,j)**6
600         app (i,j)=epp(i,j)*rri*rri 
601         bpp (i,j)=-2.0D0*epp(i,j)*rri
602         ael6(i,j)=elpp6(i,j)*4.2D0**6
603         ael3(i,j)=elpp3(i,j)*4.2D0**3
604         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
605      &                    ael6(i,j),ael3(i,j)
606         enddo
607       enddo
608 C
609 C Read side-chain interaction parameters.
610 C
611       read (isidep,*) ipot,expon
612       if (ipot.lt.1 .or. ipot.gt.5) then
613         write (iout,'(2a)') 'Error while reading SC interaction',
614      &               'potential file - unknown potential type.'
615         stop
616       endif
617       expon2=expon/2
618       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
619      & ', exponents are ',expon,2*expon 
620       goto (10,20,30,30,40) ipot
621 C----------------------- LJ potential ---------------------------------
622    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
623       if (lprint) then
624         write (iout,'(/a/)') 'Parameters of the LJ potential:'
625         write (iout,'(a/)') 'The epsilon array:'
626         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
627         write (iout,'(/a)') 'One-body parameters:'
628         write (iout,'(a,4x,a)') 'residue','sigma'
629         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
630       endif
631       goto 50
632 C----------------------- LJK potential --------------------------------
633    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
634      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
635       if (lprint) then
636         write (iout,'(/a/)') 'Parameters of the LJK potential:'
637         write (iout,'(a/)') 'The epsilon array:'
638         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
639         write (iout,'(/a)') 'One-body parameters:'
640         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
641         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
642      &        i=1,ntyp)
643       endif
644       goto 50
645 C---------------------- GB or BP potential -----------------------------
646    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
647      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
648      &  (alp(i),i=1,ntyp)
649 C For the GB potential convert sigma'**2 into chi'
650       if (ipot.eq.4) then
651         do i=1,ntyp
652           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
653         enddo
654       endif
655       if (lprint) then
656         write (iout,'(/a/)') 'Parameters of the BP potential:'
657         write (iout,'(a/)') 'The epsilon array:'
658         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
659         write (iout,'(/a)') 'One-body parameters:'
660         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
661      &       '    chip  ','    alph  '
662         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
663      &                     chip(i),alp(i),i=1,ntyp)
664       endif
665       goto 50
666 C--------------------- GBV potential -----------------------------------
667    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
668      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
669      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
670       if (lprint) then
671         write (iout,'(/a/)') 'Parameters of the GBV potential:'
672         write (iout,'(a/)') 'The epsilon array:'
673         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
674         write (iout,'(/a)') 'One-body parameters:'
675         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
676      &      's||/s_|_^2','    chip  ','    alph  '
677         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
678      &           sigii(i),chip(i),alp(i),i=1,ntyp)
679       endif
680    50 continue
681       close (isidep)
682 C-----------------------------------------------------------------------
683 C Calculate the "working" parameters of SC interactions.
684       do i=2,ntyp
685         do j=1,i-1
686           eps(i,j)=eps(j,i)
687         enddo
688       enddo
689       do i=1,ntyp
690         do j=i,ntyp
691           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
692           sigma(j,i)=sigma(i,j)
693           rs0(i,j)=dwa16*sigma(i,j)
694           rs0(j,i)=rs0(i,j)
695         enddo
696       enddo
697       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
698      & 'Working parameters of the SC interactions:',
699      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
700      & '  chi1   ','   chi2   ' 
701       do i=1,ntyp
702         do j=i,ntyp
703           epsij=eps(i,j)
704           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
705             rrij=sigma(i,j)
706           else
707             rrij=rr0(i)+rr0(j)
708           endif
709           r0(i,j)=rrij
710           r0(j,i)=rrij
711           rrij=rrij**expon
712           epsij=eps(i,j)
713           sigeps=dsign(1.0D0,epsij)
714           epsij=dabs(epsij)
715           aa(i,j)=epsij*rrij*rrij
716           bb(i,j)=-sigeps*epsij*rrij
717           aa(j,i)=aa(i,j)
718           bb(j,i)=bb(i,j)
719           if (ipot.gt.2) then
720             sigt1sq=sigma0(i)**2
721             sigt2sq=sigma0(j)**2
722             sigii1=sigii(i)
723             sigii2=sigii(j)
724             ratsig1=sigt2sq/sigt1sq
725             ratsig2=1.0D0/ratsig1
726             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
727             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
728             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
729           else
730             rsum_max=sigma(i,j)
731           endif
732 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
733             sigmaii(i,j)=rsum_max
734             sigmaii(j,i)=rsum_max 
735 c         else
736 c           sigmaii(i,j)=r0(i,j)
737 c           sigmaii(j,i)=r0(i,j)
738 c         endif
739 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
740           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
741             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
742             augm(i,j)=epsij*r_augm**(2*expon)
743 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
744             augm(j,i)=augm(i,j)
745           else
746             augm(i,j)=0.0D0
747             augm(j,i)=0.0D0
748           endif
749           if (lprint) then
750             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
751      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
752      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
753           endif
754         enddo
755       enddo
756 C
757 C Define the SC-p interaction constants
758 C
759       do i=1,20
760         do j=1,2
761           eps_scp(i,j)=-1.5d0
762           rscp(i,j)=4.0d0
763         enddo
764       enddo
765 #ifdef OLDSCP
766       do i=1,20
767 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
768 C helix formation)
769 c       aad(i,1)=0.3D0*4.0D0**12
770 C Following line for constants currently implemented
771 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
772         aad(i,1)=1.5D0*4.0D0**12
773 c       aad(i,1)=0.17D0*5.6D0**12
774         aad(i,2)=aad(i,1)
775 C "Soft" SC-p repulsion
776         bad(i,1)=0.0D0
777 C Following line for constants currently implemented
778 c       aad(i,1)=0.3D0*4.0D0**6
779 C "Hard" SC-p repulsion
780         bad(i,1)=3.0D0*4.0D0**6
781 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
782         bad(i,2)=bad(i,1)
783 c       aad(i,1)=0.0D0
784 c       aad(i,2)=0.0D0
785 c       bad(i,1)=1228.8D0
786 c       bad(i,2)=1228.8D0
787       enddo
788 #else
789 C
790 C 8/9/01 Read the SC-p interaction constants from file
791 C
792       do i=1,ntyp
793         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
794       enddo
795       do i=1,ntyp
796         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
797         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
798         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
799         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
800       enddo
801
802       if (lprint) then
803         write (iout,*) "Parameters of SC-p interactions:"
804         do i=1,20
805           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
806      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
807         enddo
808       endif
809 #endif
810 C
811 C Define the constants of the disulfide bridge
812 C
813       ebr=-5.50D0
814 c
815 c Old arbitrary potential - commented out.
816 c
817 c      dbr= 4.20D0
818 c      fbr= 3.30D0
819 c
820 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
821 c energy surface of diethyl disulfide.
822 c A. Liwo and U. Kozlowska, 11/24/03
823 c
824       D0CM = 3.78d0
825       AKCM = 15.1d0
826       AKTH = 11.0d0
827       AKCT = 12.0d0
828       V1SS =-1.08d0
829       V2SS = 7.61d0
830       V3SS = 13.7d0
831
832       write (iout,'(/a)') "Disulfide bridge parameters:"
833       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
834       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
835       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
836       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
837      & ' v3ss:',v3ss
838       return
839       end