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