increase to 20 max nset waga_homology(20)
[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 #ifdef AIX
250       call flush_(iout)
251 #else
252       call flush(iout)
253 #endif
254       endif
255 #endif
256
257 #ifdef CRYST_SC
258 C
259 C Read the parameters of the probability distribution/energy expression
260 C of the side chains.
261 C
262       do i=1,ntyp
263         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
264         if (i.eq.10) then
265           dsc_inv(i)=0.0D0
266         else
267           dsc_inv(i)=1.0D0/dsc(i)
268         endif
269         if (i.ne.10) then
270         do j=1,nlob(i)
271           do k=1,3
272             do l=1,3
273               blower(l,k,j)=0.0D0
274             enddo
275           enddo
276         enddo  
277         bsc(1,i)=0.0D0
278         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
279         do j=2,nlob(i)
280           read (irotam,*) bsc(j,i)
281           read (irotam,*) (censc(k,j,i),k=1,3),
282      &                                 ((blower(k,l,j),l=1,k),k=1,3)
283         enddo
284         do j=1,nlob(i)
285           do k=1,3
286             do l=1,k
287               akl=0.0D0
288               do m=1,3
289                 akl=akl+blower(k,m,j)*blower(l,m,j)
290               enddo
291               gaussc(k,l,j,i)=akl
292               gaussc(l,k,j,i)=akl
293             enddo
294           enddo 
295         enddo
296         endif
297       enddo
298       close (irotam)
299       if (lprint) then
300         write (iout,'(/a)') 'Parameters of side-chain local geometry'
301         do i=1,ntyp
302           nlobi=nlob(i)
303           if (nlobi.gt.0) then
304           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
305      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
306 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
307 c          write (iout,'(a,f10.4,4(16x,f10.4))')
308 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
309 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
310            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
311      &                             'log h',(bsc(j,i),j=1,nlobi)
312            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
313      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
314 c          write (iout,'(a)')
315 c         do j=1,nlobi
316 c           ind=0
317 c           do k=1,3
318 c             do l=1,k
319 c              ind=ind+1
320 c              blower(k,l,j)=gaussc(ind,j,i)
321 c             enddo
322 c           enddo
323 c         enddo
324           do k=1,3
325             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
326      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
327           enddo
328           endif
329         enddo
330       endif
331 #else
332 C
333 C Read scrot parameters for potentials determined from all-atom AM1 calculations
334 C added by Urszula Kozlowska 07/11/2007
335 C
336       do i=1,ntyp
337         read (irotam,*)
338        if (i.eq.10) then
339          read (irotam,*)
340        else
341          do j=1,65
342            read(irotam,*) sc_parmin(j,i)
343          enddo
344        endif
345       enddo
346 #endif
347       close(irotam)
348 #ifdef CRYST_TOR
349 C
350 C Read torsional parameters in old format
351 C
352       read (itorp,*) ntortyp,nterm_old
353       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
354       read (itorp,*) (itortyp(i),i=1,ntyp)
355       do i=1,ntortyp
356         do j=1,ntortyp
357           read (itorp,'(a)')
358           do k=1,nterm_old
359             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
360           enddo
361         enddo
362       enddo
363       close (itorp)
364       if (lprint) then
365         write (iout,'(/a/)') 'Torsional constants:'
366         do i=1,ntortyp
367           do j=1,ntortyp
368             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
369             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
370           enddo
371         enddo
372       endif
373 #else
374 C
375 C Read torsional parameters
376 C
377       read (itorp,*) ntortyp
378       read (itorp,*) (itortyp(i),i=1,ntyp)
379       write (iout,*) 'ntortyp',ntortyp
380       do i=1,ntortyp
381         do j=1,ntortyp
382           read (itorp,*) nterm(i,j),nlor(i,j)
383           v0ij=0.0d0
384           si=-1.0d0
385           do k=1,nterm(i,j)
386             read (itorp,*) kk,v1(k,i,j),v2(k,i,j) 
387             v0ij=v0ij+si*v1(k,i,j)
388             si=-si
389           enddo
390           do k=1,nlor(i,j)
391             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) 
392             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
393           enddo
394           v0(i,j)=v0ij
395         enddo
396       enddo
397       close (itorp)
398       if (lprint) then
399         write (iout,'(/a/)') 'Torsional constants:'
400         do i=1,ntortyp
401           do j=1,ntortyp
402             write (iout,*) 'ityp',i,' jtyp',j
403             write (iout,*) 'Fourier constants'
404             do k=1,nterm(i,j)
405               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
406             enddo
407             write (iout,*) 'Lorenz constants'
408             do k=1,nlor(i,j)
409               write (iout,'(3(1pe15.5))') 
410      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
411             enddo
412           enddo
413         enddo
414       endif
415 C
416 C 6/23/01 Read parameters for double torsionals
417 C
418       do i=1,ntortyp
419         do j=1,ntortyp
420           do k=1,ntortyp
421             read (itordp,'(3a1)') t1,t2,t3
422             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
423      &        .or. t3.ne.onelett(k)) then
424               write (iout,*) "Error in double torsional parameter file",
425      &         i,j,k,t1,t2,t3
426                stop "Error in double torsional parameter file"
427             endif
428             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
429             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
430             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
431             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
432             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
433             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
434      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
435           enddo
436         enddo
437       enddo
438       if (lprint) then
439       write (iout,*) 
440       write (iout,*) 'Constants for double torsionals'
441       do i=1,ntortyp
442         do j=1,ntortyp 
443           do k=1,ntortyp
444             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
445      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
446             write (iout,*)
447             write (iout,*) 'Single angles:'
448             do l=1,ntermd_1(i,j,k)
449               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
450      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
451      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
452             enddo
453             write (iout,*)
454             write (iout,*) 'Pairs of angles:'
455             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
456             do l=1,ntermd_2(i,j,k)
457               write (iout,'(i5,20f10.5)') 
458      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
459             enddo
460             write (iout,*)
461             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
462             do l=1,ntermd_2(i,j,k)
463               write (iout,'(i5,20f10.5)') 
464      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
465             enddo
466             write (iout,*)
467           enddo
468         enddo
469       enddo
470       endif
471 #endif
472 C Read of Side-chain backbone correlation parameters
473 C Modified 11 May 2012 by Adasko
474 CCC
475 C
476       read (isccor,*) nsccortyp
477       read (isccor,*) (isccortyp(i),i=1,ntyp)
478 c      write (iout,*) 'ntortyp',ntortyp
479       maxinter=3
480 cc maxinter is maximum interaction sites
481       do l=1,maxinter
482       do i=1,nsccortyp
483         do j=1,nsccortyp
484           read (isccor,*)
485      &nterm_sccor(i,j),nlor_sccor(i,j)
486           v0ijsccor=0.0d0
487           si=-1.0d0
488
489           do k=1,nterm_sccor(i,j)
490             read (isccor,*) kk,v1sccor(k,l,i,j)
491      &    ,v2sccor(k,l,i,j)
492             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
493             si=-si
494           enddo
495           do k=1,nlor_sccor(i,j)
496             read (isccor,*) kk,vlor1sccor(k,i,j),
497      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
498             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
499      &(1+vlor3sccor(k,i,j)**2)
500           enddo
501           v0sccor(i,j)=v0ijsccor
502         enddo
503       enddo
504       enddo
505       close (isccor)
506
507       if (lprint) then
508         write (iout,'(/a/)') 'Torsional constants:'
509         do l=1,maxinter
510         do i=1,nsccortyp
511           do j=1,nsccortyp
512             write (iout,*) 'ityp',i,' jtyp',j
513             write (iout,*) 'Fourier constants'
514             do k=1,nterm_sccor(i,j)
515               write (iout,'(2(1pe15.5))')
516      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
517             enddo
518             write (iout,*) 'Lorenz constants'
519             do k=1,nlor_sccor(i,j)
520               write (iout,'(3(1pe15.5))')
521      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
522             enddo
523           enddo
524         enddo
525         enddo
526       endif
527
528 C
529 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
530 C         interaction energy of the Gly, Ala, and Pro prototypes.
531 C
532       read (ifourier,*) nloctyp
533       do i=1,nloctyp
534         read (ifourier,*)
535         read (ifourier,*) (b(ii,i),ii=1,13)
536         if (lprint) then
537         write (iout,*) 'Type',i
538         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
539         endif
540         B1(1,i)  = b(3,i)
541         B1(2,i)  = b(5,i)
542         B1tilde(1,i) = b(3,i)
543         B1tilde(2,i) =-b(5,i) 
544         B2(1,i)  = b(2,i)
545         B2(2,i)  = b(4,i)
546         CC(1,1,i)= b(7,i)
547         CC(2,2,i)=-b(7,i)
548         CC(2,1,i)= b(9,i)
549         CC(1,2,i)= b(9,i)
550         Ctilde(1,1,i)=b(7,i)
551         Ctilde(1,2,i)=b(9,i)
552         Ctilde(2,1,i)=-b(9,i)
553         Ctilde(2,2,i)=b(7,i)
554         DD(1,1,i)= b(6,i)
555         DD(2,2,i)=-b(6,i)
556         DD(2,1,i)= b(8,i)
557         DD(1,2,i)= b(8,i)
558         Dtilde(1,1,i)=b(6,i)
559         Dtilde(1,2,i)=b(8,i)
560         Dtilde(2,1,i)=-b(8,i)
561         Dtilde(2,2,i)=b(6,i)
562         EE(1,1,i)= b(10,i)+b(11,i)
563         EE(2,2,i)=-b(10,i)+b(11,i)
564         EE(2,1,i)= b(12,i)-b(13,i)
565         EE(1,2,i)= b(12,i)+b(13,i)
566       enddo
567       if (lprint) then
568       do i=1,nloctyp
569         write (iout,*) 'Type',i
570         write (iout,*) 'B1'
571 c        write (iout,'(f10.5)') B1(:,i)
572         write(iout,*) B1(1,i),B1(2,i)
573         write (iout,*) 'B2'
574 c        write (iout,'(f10.5)') B2(:,i)
575         write(iout,*) B2(1,i),B2(2,i)
576         write (iout,*) 'CC'
577         do j=1,2
578           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
579         enddo
580         write(iout,*) 'DD'
581         do j=1,2
582           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
583         enddo
584         write(iout,*) 'EE'
585         do j=1,2
586           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
587         enddo
588       enddo
589       endif
590
591 C Read electrostatic-interaction parameters
592 C
593       if (lprint) then
594         write (iout,'(/a)') 'Electrostatic interaction constants:'
595         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
596      &            'IT','JT','APP','BPP','AEL6','AEL3'
597       endif
598       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
599       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
600       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
601       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
602       close (ielep)
603       do i=1,2
604         do j=1,2
605         rri=rpp(i,j)**6
606         app (i,j)=epp(i,j)*rri*rri 
607         bpp (i,j)=-2.0D0*epp(i,j)*rri
608         ael6(i,j)=elpp6(i,j)*4.2D0**6
609         ael3(i,j)=elpp3(i,j)*4.2D0**3
610         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
611      &                    ael6(i,j),ael3(i,j)
612         enddo
613       enddo
614 C
615 C Read side-chain interaction parameters.
616 C
617       read (isidep,*) ipot,expon
618       if (ipot.lt.1 .or. ipot.gt.5) then
619         write (iout,'(2a)') 'Error while reading SC interaction',
620      &               'potential file - unknown potential type.'
621         stop
622       endif
623       expon2=expon/2
624       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
625      & ', exponents are ',expon,2*expon 
626       goto (10,20,30,30,40) ipot
627 C----------------------- LJ potential ---------------------------------
628    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
629       if (lprint) then
630         write (iout,'(/a/)') 'Parameters of the LJ potential:'
631         write (iout,'(a/)') 'The epsilon array:'
632         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
633         write (iout,'(/a)') 'One-body parameters:'
634         write (iout,'(a,4x,a)') 'residue','sigma'
635         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
636       endif
637       goto 50
638 C----------------------- LJK potential --------------------------------
639    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
640      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
641       if (lprint) then
642         write (iout,'(/a/)') 'Parameters of the LJK potential:'
643         write (iout,'(a/)') 'The epsilon array:'
644         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
645         write (iout,'(/a)') 'One-body parameters:'
646         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
647         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
648      &        i=1,ntyp)
649       endif
650       goto 50
651 C---------------------- GB or BP potential -----------------------------
652    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
653      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
654      &  (alp(i),i=1,ntyp)
655 C For the GB potential convert sigma'**2 into chi'
656       if (ipot.eq.4) then
657         do i=1,ntyp
658           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
659         enddo
660       endif
661       if (lprint) then
662         write (iout,'(/a/)') 'Parameters of the BP potential:'
663         write (iout,'(a/)') 'The epsilon array:'
664         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
665         write (iout,'(/a)') 'One-body parameters:'
666         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
667      &       '    chip  ','    alph  '
668         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
669      &                     chip(i),alp(i),i=1,ntyp)
670       endif
671       goto 50
672 C--------------------- GBV potential -----------------------------------
673    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
674      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
675      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
676       if (lprint) then
677         write (iout,'(/a/)') 'Parameters of the GBV potential:'
678         write (iout,'(a/)') 'The epsilon array:'
679         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
680         write (iout,'(/a)') 'One-body parameters:'
681         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
682      &      's||/s_|_^2','    chip  ','    alph  '
683         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
684      &           sigii(i),chip(i),alp(i),i=1,ntyp)
685       endif
686    50 continue
687       close (isidep)
688 C-----------------------------------------------------------------------
689 C Calculate the "working" parameters of SC interactions.
690       do i=2,ntyp
691         do j=1,i-1
692           eps(i,j)=eps(j,i)
693         enddo
694       enddo
695       do i=1,ntyp
696         do j=i,ntyp
697           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
698           sigma(j,i)=sigma(i,j)
699           rs0(i,j)=dwa16*sigma(i,j)
700           rs0(j,i)=rs0(i,j)
701         enddo
702       enddo
703       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
704      & 'Working parameters of the SC interactions:',
705      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
706      & '  chi1   ','   chi2   ' 
707       do i=1,ntyp
708         do j=i,ntyp
709           epsij=eps(i,j)
710           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
711             rrij=sigma(i,j)
712           else
713             rrij=rr0(i)+rr0(j)
714           endif
715           r0(i,j)=rrij
716           r0(j,i)=rrij
717           rrij=rrij**expon
718           epsij=eps(i,j)
719           sigeps=dsign(1.0D0,epsij)
720           epsij=dabs(epsij)
721           aa(i,j)=epsij*rrij*rrij
722           bb(i,j)=-sigeps*epsij*rrij
723           aa(j,i)=aa(i,j)
724           bb(j,i)=bb(i,j)
725           if (ipot.gt.2) then
726             sigt1sq=sigma0(i)**2
727             sigt2sq=sigma0(j)**2
728             sigii1=sigii(i)
729             sigii2=sigii(j)
730             ratsig1=sigt2sq/sigt1sq
731             ratsig2=1.0D0/ratsig1
732             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
733             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
734             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
735           else
736             rsum_max=sigma(i,j)
737           endif
738 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
739             sigmaii(i,j)=rsum_max
740             sigmaii(j,i)=rsum_max 
741 c         else
742 c           sigmaii(i,j)=r0(i,j)
743 c           sigmaii(j,i)=r0(i,j)
744 c         endif
745 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
746           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
747             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
748             augm(i,j)=epsij*r_augm**(2*expon)
749 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
750             augm(j,i)=augm(i,j)
751           else
752             augm(i,j)=0.0D0
753             augm(j,i)=0.0D0
754           endif
755           if (lprint) then
756             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
757      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
758      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
759           endif
760         enddo
761       enddo
762 C
763 C Define the SC-p interaction constants
764 C
765       do i=1,20
766         do j=1,2
767           eps_scp(i,j)=-1.5d0
768           rscp(i,j)=4.0d0
769         enddo
770       enddo
771 #ifdef OLDSCP
772       do i=1,20
773 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
774 C helix formation)
775 c       aad(i,1)=0.3D0*4.0D0**12
776 C Following line for constants currently implemented
777 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
778         aad(i,1)=1.5D0*4.0D0**12
779 c       aad(i,1)=0.17D0*5.6D0**12
780         aad(i,2)=aad(i,1)
781 C "Soft" SC-p repulsion
782         bad(i,1)=0.0D0
783 C Following line for constants currently implemented
784 c       aad(i,1)=0.3D0*4.0D0**6
785 C "Hard" SC-p repulsion
786         bad(i,1)=3.0D0*4.0D0**6
787 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
788         bad(i,2)=bad(i,1)
789 c       aad(i,1)=0.0D0
790 c       aad(i,2)=0.0D0
791 c       bad(i,1)=1228.8D0
792 c       bad(i,2)=1228.8D0
793       enddo
794 #else
795 C
796 C 8/9/01 Read the SC-p interaction constants from file
797 C
798       do i=1,ntyp
799         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
800       enddo
801       do i=1,ntyp
802         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
803         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
804         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
805         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
806       enddo
807
808       if (lprint) then
809         write (iout,*) "Parameters of SC-p interactions:"
810         do i=1,20
811           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
812      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
813         enddo
814       endif
815 #endif
816 C
817 C Define the constants of the disulfide bridge
818 C
819 C      ebr=-5.50D0
820 c
821 c Old arbitrary potential - commented out.
822 c
823 c      dbr= 4.20D0
824 c      fbr= 3.30D0
825 c
826 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
827 c energy surface of diethyl disulfide.
828 c A. Liwo and U. Kozlowska, 11/24/03
829 c
830 C      D0CM = 3.78d0
831 C      AKCM = 15.1d0
832 C      AKTH = 11.0d0
833 C      AKCT = 12.0d0
834 C      V1SS =-1.08d0
835 C      V2SS = 7.61d0
836 C      V3SS = 13.7d0
837       return
838       end