arcos = 0.5D0*(PI-DSIGN(1.0D0,X)*PI)
[unres.git] / source / cluster / wham / src-M / 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
469 C Read of Side-chain backbone correlation parameters
470 C Modified 11 May 2012 by Adasko
471 CCC
472       read (isccor,*) nsccortyp
473       read (isccor,*) (isccortyp(i),i=1,ntyp)
474       do i=-ntyp,-1
475         isccortyp(i)=-isccortyp(-i)
476       enddo
477       iscprol=isccortyp(20)
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           v0ijsccor1=0.0d0
488           v0ijsccor2=0.0d0
489           v0ijsccor3=0.0d0
490           si=-1.0d0
491           nterm_sccor(-i,j)=nterm_sccor(i,j)
492           nterm_sccor(-i,-j)=nterm_sccor(i,j)
493           nterm_sccor(i,-j)=nterm_sccor(i,j)
494           do k=1,nterm_sccor(i,j)
495             read (isccor,*) kk,v1sccor(k,l,i,j)
496      &    ,v2sccor(k,l,i,j)
497             if (j.eq.iscprol) then
498              if (i.eq.isccortyp(10)) then
499              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
500              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
501              else
502              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
503      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
504              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
505      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
506              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
507              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
508              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
509              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
510              endif
511             else
512              if (i.eq.isccortyp(10)) then
513              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
514              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
515              else
516                if (j.eq.isccortyp(10)) then
517              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
518              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
519                else
520              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
521              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
522              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
523              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
524              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
525              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
526                 endif
527                endif
528             endif
529             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
530             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
531             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
532             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
533             si=-si
534            enddo
535           do k=1,nlor_sccor(i,j)
536             read (isccor,*) kk,vlor1sccor(k,i,j),
537      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
538             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
539      &(1+vlor3sccor(k,i,j)**2)
540           enddo
541           v0sccor(l,i,j)=v0ijsccor
542           v0sccor(l,-i,j)=v0ijsccor1
543           v0sccor(l,i,-j)=v0ijsccor2
544           v0sccor(l,-i,-j)=v0ijsccor3             
545          enddo
546         enddo
547       enddo
548       close (isccor)
549       if (lprint) then
550         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
551         do i=1,nsccortyp
552           do j=1,nsccortyp
553             write (iout,*) 'ityp',i,' jtyp',j
554             write (iout,*) 'Fourier constants'
555             do k=1,nterm_sccor(i,j)
556               write (iout,'(2(1pe15.5))')
557      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
558             enddo
559             write (iout,*) 'Lorenz constants'
560             do k=1,nlor_sccor(i,j)
561               write (iout,'(3(1pe15.5))')
562      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
563             enddo 
564           enddo
565         enddo
566       endif
567 C
568 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
569 C         interaction energy of the Gly, Ala, and Pro prototypes.
570 C
571       read (ifourier,*) nloctyp
572       do i=1,nloctyp
573         read (ifourier,*)
574         read (ifourier,*) (b(ii,i),ii=1,13)
575         if (lprint) then
576         write (iout,*) 'Type',i
577         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
578         endif
579         B1(1,i)  = b(3,i)
580         B1(2,i)  = b(5,i)
581         B1tilde(1,i) = b(3,i)
582         B1tilde(2,i) =-b(5,i) 
583         B2(1,i)  = b(2,i)
584         B2(2,i)  = b(4,i)
585         CC(1,1,i)= b(7,i)
586         CC(2,2,i)=-b(7,i)
587         CC(2,1,i)= b(9,i)
588         CC(1,2,i)= b(9,i)
589         Ctilde(1,1,i)=b(7,i)
590         Ctilde(1,2,i)=b(9,i)
591         Ctilde(2,1,i)=-b(9,i)
592         Ctilde(2,2,i)=b(7,i)
593         DD(1,1,i)= b(6,i)
594         DD(2,2,i)=-b(6,i)
595         DD(2,1,i)= b(8,i)
596         DD(1,2,i)= b(8,i)
597         Dtilde(1,1,i)=b(6,i)
598         Dtilde(1,2,i)=b(8,i)
599         Dtilde(2,1,i)=-b(8,i)
600         Dtilde(2,2,i)=b(6,i)
601         EE(1,1,i)= b(10,i)+b(11,i)
602         EE(2,2,i)=-b(10,i)+b(11,i)
603         EE(2,1,i)= b(12,i)-b(13,i)
604         EE(1,2,i)= b(12,i)+b(13,i)
605       enddo
606       if (lprint) then
607       do i=1,nloctyp
608         write (iout,*) 'Type',i
609         write (iout,*) 'B1'
610 c        write (iout,'(f10.5)') B1(:,i)
611         write(iout,*) B1(1,i),B1(2,i)
612         write (iout,*) 'B2'
613 c        write (iout,'(f10.5)') B2(:,i)
614         write(iout,*) B2(1,i),B2(2,i)
615         write (iout,*) 'CC'
616         do j=1,2
617           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
618         enddo
619         write(iout,*) 'DD'
620         do j=1,2
621           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
622         enddo
623         write(iout,*) 'EE'
624         do j=1,2
625           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
626         enddo
627       enddo
628       endif
629
630 C Read electrostatic-interaction parameters
631 C
632       if (lprint) then
633         write (iout,'(/a)') 'Electrostatic interaction constants:'
634         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
635      &            'IT','JT','APP','BPP','AEL6','AEL3'
636       endif
637       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
638       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
639       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
640       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
641       close (ielep)
642       do i=1,2
643         do j=1,2
644         rri=rpp(i,j)**6
645         app (i,j)=epp(i,j)*rri*rri 
646         bpp (i,j)=-2.0D0*epp(i,j)*rri
647         ael6(i,j)=elpp6(i,j)*4.2D0**6
648         ael3(i,j)=elpp3(i,j)*4.2D0**3
649         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
650      &                    ael6(i,j),ael3(i,j)
651         enddo
652       enddo
653 C
654 C Read side-chain interaction parameters.
655 C
656       read (isidep,*) ipot,expon
657       if (ipot.lt.1 .or. ipot.gt.5) then
658         write (iout,'(2a)') 'Error while reading SC interaction',
659      &               'potential file - unknown potential type.'
660         stop
661       endif
662       expon2=expon/2
663       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
664      & ', exponents are ',expon,2*expon 
665       goto (10,20,30,30,40) ipot
666 C----------------------- LJ potential ---------------------------------
667    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
668       if (lprint) then
669         write (iout,'(/a/)') 'Parameters of the LJ potential:'
670         write (iout,'(a/)') 'The epsilon array:'
671         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
672         write (iout,'(/a)') 'One-body parameters:'
673         write (iout,'(a,4x,a)') 'residue','sigma'
674         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
675       endif
676       goto 50
677 C----------------------- LJK potential --------------------------------
678    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
679      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
680       if (lprint) then
681         write (iout,'(/a/)') 'Parameters of the LJK potential:'
682         write (iout,'(a/)') 'The epsilon array:'
683         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
684         write (iout,'(/a)') 'One-body parameters:'
685         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
686         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
687      &        i=1,ntyp)
688       endif
689       goto 50
690 C---------------------- GB or BP potential -----------------------------
691    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
692      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
693      &  (alp(i),i=1,ntyp)
694 C For the GB potential convert sigma'**2 into chi'
695       if (ipot.eq.4) then
696         do i=1,ntyp
697           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
698         enddo
699       endif
700       if (lprint) then
701         write (iout,'(/a/)') 'Parameters of the BP potential:'
702         write (iout,'(a/)') 'The epsilon array:'
703         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
704         write (iout,'(/a)') 'One-body parameters:'
705         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
706      &       '    chip  ','    alph  '
707         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
708      &                     chip(i),alp(i),i=1,ntyp)
709       endif
710       goto 50
711 C--------------------- GBV potential -----------------------------------
712    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
713      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
714      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
715       if (lprint) then
716         write (iout,'(/a/)') 'Parameters of the GBV potential:'
717         write (iout,'(a/)') 'The epsilon array:'
718         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
719         write (iout,'(/a)') 'One-body parameters:'
720         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
721      &      's||/s_|_^2','    chip  ','    alph  '
722         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
723      &           sigii(i),chip(i),alp(i),i=1,ntyp)
724       endif
725    50 continue
726       close (isidep)
727 C-----------------------------------------------------------------------
728 C Calculate the "working" parameters of SC interactions.
729       do i=2,ntyp
730         do j=1,i-1
731           eps(i,j)=eps(j,i)
732         enddo
733       enddo
734       do i=1,ntyp
735         do j=i,ntyp
736           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
737           sigma(j,i)=sigma(i,j)
738           rs0(i,j)=dwa16*sigma(i,j)
739           rs0(j,i)=rs0(i,j)
740         enddo
741       enddo
742       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
743      & 'Working parameters of the SC interactions:',
744      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
745      & '  chi1   ','   chi2   ' 
746       do i=1,ntyp
747         do j=i,ntyp
748           epsij=eps(i,j)
749           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
750             rrij=sigma(i,j)
751           else
752             rrij=rr0(i)+rr0(j)
753           endif
754           r0(i,j)=rrij
755           r0(j,i)=rrij
756           rrij=rrij**expon
757           epsij=eps(i,j)
758           sigeps=dsign(1.0D0,epsij)
759           epsij=dabs(epsij)
760           aa(i,j)=epsij*rrij*rrij
761           bb(i,j)=-sigeps*epsij*rrij
762           aa(j,i)=aa(i,j)
763           bb(j,i)=bb(i,j)
764           if (ipot.gt.2) then
765             sigt1sq=sigma0(i)**2
766             sigt2sq=sigma0(j)**2
767             sigii1=sigii(i)
768             sigii2=sigii(j)
769             ratsig1=sigt2sq/sigt1sq
770             ratsig2=1.0D0/ratsig1
771             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
772             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
773             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
774           else
775             rsum_max=sigma(i,j)
776           endif
777 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
778             sigmaii(i,j)=rsum_max
779             sigmaii(j,i)=rsum_max 
780 c         else
781 c           sigmaii(i,j)=r0(i,j)
782 c           sigmaii(j,i)=r0(i,j)
783 c         endif
784 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
785           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
786             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
787             augm(i,j)=epsij*r_augm**(2*expon)
788 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
789             augm(j,i)=augm(i,j)
790           else
791             augm(i,j)=0.0D0
792             augm(j,i)=0.0D0
793           endif
794           if (lprint) then
795             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
796      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
797      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
798           endif
799         enddo
800       enddo
801 C
802 C Define the SC-p interaction constants
803 C
804       do i=1,20
805         do j=1,2
806           eps_scp(i,j)=-1.5d0
807           rscp(i,j)=4.0d0
808         enddo
809       enddo
810 #ifdef OLDSCP
811       do i=1,20
812 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
813 C helix formation)
814 c       aad(i,1)=0.3D0*4.0D0**12
815 C Following line for constants currently implemented
816 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
817         aad(i,1)=1.5D0*4.0D0**12
818 c       aad(i,1)=0.17D0*5.6D0**12
819         aad(i,2)=aad(i,1)
820 C "Soft" SC-p repulsion
821         bad(i,1)=0.0D0
822 C Following line for constants currently implemented
823 c       aad(i,1)=0.3D0*4.0D0**6
824 C "Hard" SC-p repulsion
825         bad(i,1)=3.0D0*4.0D0**6
826 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
827         bad(i,2)=bad(i,1)
828 c       aad(i,1)=0.0D0
829 c       aad(i,2)=0.0D0
830 c       bad(i,1)=1228.8D0
831 c       bad(i,2)=1228.8D0
832       enddo
833 #else
834 C
835 C 8/9/01 Read the SC-p interaction constants from file
836 C
837       do i=1,ntyp
838         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
839       enddo
840       do i=1,ntyp
841         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
842         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
843         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
844         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
845       enddo
846
847       if (lprint) then
848         write (iout,*) "Parameters of SC-p interactions:"
849         do i=1,20
850           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
851      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
852         enddo
853       endif
854 #endif
855 C
856 C Define the constants of the disulfide bridge
857 C
858       ebr=-5.50D0
859 c
860 c Old arbitrary potential - commented out.
861 c
862 c      dbr= 4.20D0
863 c      fbr= 3.30D0
864 c
865 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
866 c energy surface of diethyl disulfide.
867 c A. Liwo and U. Kozlowska, 11/24/03
868 c
869       D0CM = 3.78d0
870       AKCM = 15.1d0
871       AKTH = 11.0d0
872       AKCT = 12.0d0
873       V1SS =-1.08d0
874       V2SS = 7.61d0
875       V3SS = 13.7d0
876
877       write (iout,'(/a)') "Disulfide bridge parameters:"
878       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
879       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
880       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
881       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
882      & ' v3ss:',v3ss
883       return
884       end