8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
874 if (oldion.eq.1) then
876 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
877 print *,msc(i,5),restok(i,5)
881 read (iion,*) ncatprotparm
882 allocate(catprm(ncatprotparm,4))
884 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
888 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
892 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
895 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
896 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
900 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
901 restyp(i,5),"-",restyp(j,2)
904 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
905 !----------------------------------------------------
906 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
907 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
908 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
909 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
910 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
911 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
921 allocate(liptranene(ntyp))
922 !C reading lipid parameters
923 write (iout,*) "iliptranpar",iliptranpar
925 read(iliptranpar,*) pepliptran
928 read(iliptranpar,*) liptranene(i)
929 print *,liptranene(i)
935 ! Read the parameters of the probability distribution/energy expression
936 ! of the virtual-bond valence angles theta
939 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
940 (bthet(j,i,1,1),j=1,2)
941 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
942 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
943 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
947 athet(1,i,1,-1)=athet(1,i,1,1)
948 athet(2,i,1,-1)=athet(2,i,1,1)
949 bthet(1,i,1,-1)=-bthet(1,i,1,1)
950 bthet(2,i,1,-1)=-bthet(2,i,1,1)
951 athet(1,i,-1,1)=-athet(1,i,1,1)
952 athet(2,i,-1,1)=-athet(2,i,1,1)
953 bthet(1,i,-1,1)=bthet(1,i,1,1)
954 bthet(2,i,-1,1)=bthet(2,i,1,1)
958 athet(1,i,-1,-1)=athet(1,-i,1,1)
959 athet(2,i,-1,-1)=-athet(2,-i,1,1)
960 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
961 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
962 athet(1,i,-1,1)=athet(1,-i,1,1)
963 athet(2,i,-1,1)=-athet(2,-i,1,1)
964 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
965 bthet(2,i,-1,1)=bthet(2,-i,1,1)
966 athet(1,i,1,-1)=-athet(1,-i,1,1)
967 athet(2,i,1,-1)=athet(2,-i,1,1)
968 bthet(1,i,1,-1)=bthet(1,-i,1,1)
969 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
974 polthet(j,i)=polthet(j,-i)
977 gthet(j,i)=gthet(j,-i)
985 'Parameters of the virtual-bond valence angles:'
986 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
987 ' ATHETA0 ',' A1 ',' A2 ',&
990 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
991 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
993 write (iout,'(/a/9x,5a/79(1h-))') &
994 'Parameters of the expression for sigma(theta_c):',&
995 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
996 ' ALPH3 ',' SIGMA0C '
998 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
999 (polthet(j,i),j=0,3),sigc0(i)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the second gaussian:',&
1003 ' THETA0 ',' SIGMA0 ',' G1 ',&
1006 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1007 sig0(i),(gthet(j,i),j=1,3)
1010 write (iout,'(a)') &
1011 'Parameters of the virtual-bond valence angles:'
1012 write (iout,'(/a/9x,5a/79(1h-))') &
1013 'Coefficients of expansion',&
1014 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1015 ' b1*10^1 ',' b2*10^1 '
1017 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1018 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1019 (10*bthet(j,i,1,1),j=1,2)
1021 write (iout,'(/a/9x,5a/79(1h-))') &
1022 'Parameters of the expression for sigma(theta_c):',&
1023 ' alpha0 ',' alph1 ',' alph2 ',&
1024 ' alhp3 ',' sigma0c '
1026 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1027 (polthet(j,i),j=0,3),sigc0(i)
1029 write (iout,'(/a/9x,5a/79(1h-))') &
1030 'Parameters of the second gaussian:',&
1031 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1034 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1035 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1041 ! Read the parameters of Utheta determined from ab initio surfaces
1042 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1044 IF (tor_mode.eq.0) THEN
1045 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1046 ntheterm3,nsingle,ndouble
1047 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1049 !----------------------------------------------------
1050 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1051 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1052 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1053 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1054 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1055 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1056 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1057 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1058 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1059 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1060 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1061 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1062 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1063 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1064 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1065 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1066 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1067 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1068 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1069 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1070 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1071 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1072 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1073 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1075 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1077 ithetyp(i)=-ithetyp(-i)
1080 aa0thet(:,:,:,:)=0.0d0
1081 aathet(:,:,:,:,:)=0.0d0
1082 bbthet(:,:,:,:,:,:)=0.0d0
1083 ccthet(:,:,:,:,:,:)=0.0d0
1084 ddthet(:,:,:,:,:,:)=0.0d0
1085 eethet(:,:,:,:,:,:)=0.0d0
1086 ffthet(:,:,:,:,:,:,:)=0.0d0
1087 ggthet(:,:,:,:,:,:,:)=0.0d0
1089 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1091 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1092 ! VAR:1=non-glicyne non-proline 2=proline
1093 ! VAR:negative values for D-aminoacid
1095 do j=-nthetyp,nthetyp
1096 do k=-nthetyp,nthetyp
1097 read (ithep,'(6a)',end=111,err=111) res1
1098 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1099 ! VAR: aa0thet is variable describing the average value of Foureir
1100 ! VAR: expansion series
1101 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1102 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1103 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1104 read (ithep,*,end=111,err=111) &
1105 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1106 read (ithep,*,end=111,err=111) &
1107 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1108 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1109 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1110 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1112 read (ithep,*,end=111,err=111) &
1113 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1114 ffthet(lll,llll,ll,i,j,k,iblock),&
1115 ggthet(llll,lll,ll,i,j,k,iblock),&
1116 ggthet(lll,llll,ll,i,j,k,iblock),&
1117 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1122 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1123 ! coefficients of theta-and-gamma-dependent terms are zero.
1124 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1125 ! RECOMENTDED AFTER VERSION 3.3)
1129 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1130 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1132 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1133 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1136 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1138 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1141 ! AND COMMENT THE LOOPS BELOW
1145 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1146 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1148 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1149 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1152 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1154 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1159 ! Substitution for D aminoacids from symmetry.
1162 do j=-nthetyp,nthetyp
1163 do k=-nthetyp,nthetyp
1164 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1166 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1170 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1171 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1172 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1173 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1179 ffthet(llll,lll,ll,i,j,k,iblock)= &
1180 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1181 ffthet(lll,llll,ll,i,j,k,iblock)= &
1182 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1183 ggthet(llll,lll,ll,i,j,k,iblock)= &
1184 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1185 ggthet(lll,llll,ll,i,j,k,iblock)= &
1186 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1195 ! Control printout of the coefficients of virtual-bond-angle potentials
1198 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1203 write (iout,'(//4a)') &
1204 'Type ',onelett(i),onelett(j),onelett(k)
1205 write (iout,'(//a,10x,a)') " l","a[l]"
1206 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1207 write (iout,'(i2,1pe15.5)') &
1208 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1210 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1211 "b",l,"c",l,"d",l,"e",l
1213 write (iout,'(i2,4(1pe15.5))') m,&
1214 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1215 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1219 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1220 "f+",l,"f-",l,"g+",l,"g-",l
1223 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1224 ffthet(n,m,l,i,j,k,iblock),&
1225 ffthet(m,n,l,i,j,k,iblock),&
1226 ggthet(n,m,l,i,j,k,iblock),&
1227 ggthet(m,n,l,i,j,k,iblock)
1238 !C here will be the apropriate recalibrating for D-aminoacid
1239 read (ithep,*,end=121,err=121) nthetyp
1240 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1241 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1242 do i=-nthetyp+1,nthetyp-1
1243 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1244 do j=0,nbend_kcc_Tb(i)
1245 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1249 write (iout,'(a)') &
1250 "Parameters of the valence-only potentials"
1251 do i=-nthetyp+1,nthetyp-1
1252 write (iout,'(2a)') "Type ",toronelet(i)
1253 do k=0,nbend_kcc_Tb(i)
1254 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1260 write (2,*) "Start reading THETA_PDB",ithep_pdb
1262 ! write (2,*) 'i=',i
1263 read (ithep_pdb,*,err=111,end=111) &
1264 a0thet(i),(athet(j,i,1,1),j=1,2),&
1265 (bthet(j,i,1,1),j=1,2)
1266 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1267 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1268 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1269 sigc0(i)=sigc0(i)**2
1272 athet(1,i,1,-1)=athet(1,i,1,1)
1273 athet(2,i,1,-1)=athet(2,i,1,1)
1274 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1275 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1276 athet(1,i,-1,1)=-athet(1,i,1,1)
1277 athet(2,i,-1,1)=-athet(2,i,1,1)
1278 bthet(1,i,-1,1)=bthet(1,i,1,1)
1279 bthet(2,i,-1,1)=bthet(2,i,1,1)
1282 a0thet(i)=a0thet(-i)
1283 athet(1,i,-1,-1)=athet(1,-i,1,1)
1284 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1285 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1286 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1287 athet(1,i,-1,1)=athet(1,-i,1,1)
1288 athet(2,i,-1,1)=-athet(2,-i,1,1)
1289 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1290 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1291 athet(1,i,1,-1)=-athet(1,-i,1,1)
1292 athet(2,i,1,-1)=athet(2,-i,1,1)
1293 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1294 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1295 theta0(i)=theta0(-i)
1299 polthet(j,i)=polthet(j,-i)
1302 gthet(j,i)=gthet(j,-i)
1305 write (2,*) "End reading THETA_PDB"
1309 !--------------- Reading theta parameters for nucleic acid-------
1310 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1311 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1312 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1313 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1314 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1315 nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1317 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1318 nthetyp_nucl+1,nthetyp_nucl+1))
1319 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1320 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1321 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1322 nthetyp_nucl+1,nthetyp_nucl+1))
1323 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1324 nthetyp_nucl+1,nthetyp_nucl+1))
1325 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1326 nthetyp_nucl+1,nthetyp_nucl+1))
1327 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1328 nthetyp_nucl+1,nthetyp_nucl+1))
1329 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1330 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1331 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1332 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1333 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1334 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1336 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1337 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1339 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1341 aa0thet_nucl(:,:,:)=0.0d0
1342 aathet_nucl(:,:,:,:)=0.0d0
1343 bbthet_nucl(:,:,:,:,:)=0.0d0
1344 ccthet_nucl(:,:,:,:,:)=0.0d0
1345 ddthet_nucl(:,:,:,:,:)=0.0d0
1346 eethet_nucl(:,:,:,:,:)=0.0d0
1347 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1348 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1353 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1354 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1355 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1356 read (ithep_nucl,*,end=111,err=111) &
1357 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1358 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1359 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1360 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1361 read (ithep_nucl,*,end=111,err=111) &
1362 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1363 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1364 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1369 !-------------------------------------------
1370 allocate(nlob(ntyp1)) !(ntyp1)
1371 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1372 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1373 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1380 gaussc(:,:,:,:)=0.0D0
1384 ! Read the parameters of the probability distribution/energy expression
1385 ! of the side chains.
1388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1392 dsc_inv(i)=1.0D0/dsc(i)
1403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1404 ((blower(k,l,1),l=1,k),k=1,3)
1405 censc(1,1,-i)=censc(1,1,i)
1406 censc(2,1,-i)=censc(2,1,i)
1407 censc(3,1,-i)=-censc(3,1,i)
1409 read (irotam,*,end=112,err=112) bsc(j,i)
1410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1411 ((blower(k,l,j),l=1,k),k=1,3)
1412 censc(1,j,-i)=censc(1,j,i)
1413 censc(2,j,-i)=censc(2,j,i)
1414 censc(3,j,-i)=-censc(3,j,i)
1415 ! BSC is amplitude of Gaussian
1422 akl=akl+blower(k,m,j)*blower(l,m,j)
1426 if (((k.eq.3).and.(l.ne.3)) &
1427 .or.((l.eq.3).and.(k.ne.3))) then
1428 gaussc(k,l,j,-i)=-akl
1429 gaussc(l,k,j,-i)=-akl
1431 gaussc(k,l,j,-i)=akl
1432 gaussc(l,k,j,-i)=akl
1441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1444 if (nlobi.gt.0) then
1446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1447 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1449 'log h',(bsc(j,i),j=1,nlobi)
1450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1451 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1454 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1458 write (iout,'(a,f10.4,4(16x,f10.4))') &
1459 'Center ',(bsc(j,i),j=1,nlobi)
1460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1469 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1470 ! added by Urszula Kozlowska 07/11/2007
1472 !el Maximum number of SC local term fitting function coefficiants
1473 !el integer,parameter :: maxsccoef=65
1475 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1478 read (irotam,*,end=112,err=112)
1480 read (irotam,*,end=112,err=112)
1483 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1487 !---------reading nucleic acid parameters for rotamers-------------------
1488 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1489 do i=1,ntyp_molec(2)
1490 read (irotam_nucl,*,end=112,err=112)
1492 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1498 write (iout,*) "Base rotamer parameters"
1499 do i=1,ntyp_molec(2)
1500 write (iout,'(a)') restyp(i,2)
1501 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1506 ! Read the parameters of the probability distribution/energy expression
1507 ! of the side chains.
1509 write (2,*) "Start reading ROTAM_PDB"
1511 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1515 dsc_inv(i)=1.0D0/dsc(i)
1526 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1527 ((blower(k,l,1),l=1,k),k=1,3)
1529 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1530 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1531 ((blower(k,l,j),l=1,k),k=1,3)
1538 akl=akl+blower(k,m,j)*blower(l,m,j)
1548 write (2,*) "End reading ROTAM_PDB"
1554 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1555 !C interaction energy of the Gly, Ala, and Pro prototypes.
1557 read (ifourier,*) nloctyp
1558 SPLIT_FOURIERTOR = nloctyp.lt.0
1559 nloctyp = iabs(nloctyp)
1560 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1561 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1562 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1563 !C allocate(ctilde(2,2,nres))
1564 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1565 !C allocate(gtb1(2,nres))
1566 !C allocate(gtb2(2,nres))
1567 !C allocate(cc(2,2,nres))
1568 !C allocate(dd(2,2,nres))
1569 !C allocate(ee(2,2,nres))
1570 !C allocate(gtcc(2,2,nres))
1571 !C allocate(gtdd(2,2,nres))
1572 !C allocate(gtee(2,2,nres))
1575 allocate(itype2loc(-ntyp1:ntyp1))
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(bnew1(3,2,-nloctyp:nloctyp))
1578 allocate(bnew2(3,2,-nloctyp:nloctyp))
1579 allocate(ccnew(3,2,-nloctyp:nloctyp))
1580 allocate(ddnew(3,2,-nloctyp:nloctyp))
1581 allocate(e0new(3,-nloctyp:nloctyp))
1582 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1583 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1584 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1585 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1586 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1587 allocate(e0newtor(3,-nloctyp:nloctyp))
1588 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1590 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1591 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1592 itype2loc(ntyp1)=nloctyp
1593 iloctyp(nloctyp)=ntyp1
1595 itype2loc(-i)=-itype2loc(i)
1598 allocate(iloctyp(-nloctyp:nloctyp))
1599 allocate(itype2loc(-ntyp1:ntyp1))
1606 iloctyp(-i)=-iloctyp(i)
1608 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1609 !c write (iout,*) "nloctyp",nloctyp,
1610 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1611 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1612 !c write (iout,*) "nloctyp",nloctyp,
1613 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1616 !c write (iout,*) "NEWCORR",i
1617 read (ifourier,*,end=115,err=115)
1620 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1623 !c write (iout,*) "NEWCORR BNEW1"
1624 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1630 !c write (iout,*) "NEWCORR BNEW2"
1631 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1633 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1634 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1636 !c write (iout,*) "NEWCORR CCNEW"
1637 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1639 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1640 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1642 !c write (iout,*) "NEWCORR DDNEW"
1643 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1647 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1651 !c write (iout,*) "NEWCORR EENEW1"
1652 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1654 read (ifourier,*,end=115,err=115) e0new(ii,i)
1656 !c write (iout,*) (e0new(ii,i),ii=1,3)
1658 !c write (iout,*) "NEWCORR EENEW"
1661 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1662 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1663 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1664 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1669 bnew1(ii,1,-i)= bnew1(ii,1,i)
1670 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1671 bnew2(ii,1,-i)= bnew2(ii,1,i)
1672 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1675 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1676 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1677 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1678 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1679 ccnew(ii,1,-i)=ccnew(ii,1,i)
1680 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1681 ddnew(ii,1,-i)=ddnew(ii,1,i)
1682 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1684 e0new(1,-i)= e0new(1,i)
1685 e0new(2,-i)=-e0new(2,i)
1686 e0new(3,-i)=-e0new(3,i)
1688 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1689 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1690 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1691 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1695 write (iout,'(a)') "Coefficients of the multibody terms"
1696 do i=-nloctyp+1,nloctyp-1
1697 write (iout,*) "Type: ",onelet(iloctyp(i))
1698 write (iout,*) "Coefficients of the expansion of B1"
1700 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1702 write (iout,*) "Coefficients of the expansion of B2"
1704 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1706 write (iout,*) "Coefficients of the expansion of C"
1707 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1708 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1709 write (iout,*) "Coefficients of the expansion of D"
1710 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1711 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1712 write (iout,*) "Coefficients of the expansion of E"
1713 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1716 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1721 IF (SPLIT_FOURIERTOR) THEN
1723 !c write (iout,*) "NEWCORR TOR",i
1724 read (ifourier,*,end=115,err=115)
1727 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1730 !c write (iout,*) "NEWCORR BNEW1 TOR"
1731 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1737 !c write (iout,*) "NEWCORR BNEW2 TOR"
1738 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1740 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1741 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1743 !c write (iout,*) "NEWCORR CCNEW TOR"
1744 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1746 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1747 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1749 !c write (iout,*) "NEWCORR DDNEW TOR"
1750 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1754 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1758 !c write (iout,*) "NEWCORR EENEW1 TOR"
1759 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1761 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1763 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1765 !c write (iout,*) "NEWCORR EENEW TOR"
1768 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1769 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1770 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1771 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1776 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1777 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1778 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1779 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1782 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1783 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1784 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1785 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1787 e0newtor(1,-i)= e0newtor(1,i)
1788 e0newtor(2,-i)=-e0newtor(2,i)
1789 e0newtor(3,-i)=-e0newtor(3,i)
1791 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1792 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1793 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1794 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1798 write (iout,'(a)') &
1799 "Single-body coefficients of the torsional potentials"
1800 do i=-nloctyp+1,nloctyp-1
1801 write (iout,*) "Type: ",onelet(iloctyp(i))
1802 write (iout,*) "Coefficients of the expansion of B1tor"
1804 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1805 j,(bnew1tor(k,j,i),k=1,3)
1807 write (iout,*) "Coefficients of the expansion of B2tor"
1809 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1810 j,(bnew2tor(k,j,i),k=1,3)
1812 write (iout,*) "Coefficients of the expansion of Ctor"
1813 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1814 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1815 write (iout,*) "Coefficients of the expansion of Dtor"
1816 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1817 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1818 write (iout,*) "Coefficients of the expansion of Etor"
1819 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1822 write (iout,'(1hE,2i1,2f10.5)') &
1823 j,k,(eenewtor(l,j,k,i),l=1,2)
1829 do i=-nloctyp+1,nloctyp-1
1832 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1837 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1841 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1842 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1843 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1844 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1847 ENDIF !SPLIT_FOURIER_TOR
1849 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1850 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1851 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1852 allocate(b(13,-nloctyp-1:nloctyp+1))
1854 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1856 read (ifourier,*,end=115,err=115)
1857 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1859 write (iout,*) 'Type ',onelet(iloctyp(i))
1860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1874 !c B1tilde(1,i) = b(3)
1875 !c! B1tilde(2,i) =-b(5)
1876 !c! B1tilde(1,-i) =-b(3)
1877 !c! B1tilde(2,-i) =b(5)
1878 !c! b1tilde(1,i)=0.0d0
1879 !c b1tilde(2,i)=0.0d0
1884 !cc B1tilde(1,i) = b(3,i)
1885 !cc B1tilde(2,i) =-b(5,i)
1886 !c B1tilde(1,-i) =-b(3,i)
1887 !c B1tilde(2,-i) =b(5,i)
1888 !cc b1tilde(1,i)=0.0d0
1889 !cc b1tilde(2,i)=0.0d0
1890 !cc B2(1,i) = b(2,i)
1891 !cc B2(2,i) = b(4,i)
1893 !c B2(2,-i) =-b(4,i)
1897 CCold(1,1,i)= b(7,i)
1898 CCold(2,2,i)=-b(7,i)
1899 CCold(2,1,i)= b(9,i)
1900 CCold(1,2,i)= b(9,i)
1901 CCold(1,1,-i)= b(7,i)
1902 CCold(2,2,-i)=-b(7,i)
1903 CCold(2,1,-i)=-b(9,i)
1904 CCold(1,2,-i)=-b(9,i)
1909 !c Ctilde(1,1,i)= CCold(1,1,i)
1910 !c Ctilde(1,2,i)= CCold(1,2,i)
1911 !c Ctilde(2,1,i)=-CCold(2,1,i)
1912 !c Ctilde(2,2,i)=-CCold(2,2,i)
1917 !c Ctilde(1,1,i)= CCold(1,1,i)
1918 !c Ctilde(1,2,i)= CCold(1,2,i)
1919 !c Ctilde(2,1,i)=-CCold(2,1,i)
1920 !c Ctilde(2,2,i)=-CCold(2,2,i)
1922 !c Ctilde(1,1,i)=0.0d0
1923 !c Ctilde(1,2,i)=0.0d0
1924 !c Ctilde(2,1,i)=0.0d0
1925 !c Ctilde(2,2,i)=0.0d0
1926 DDold(1,1,i)= b(6,i)
1927 DDold(2,2,i)=-b(6,i)
1928 DDold(2,1,i)= b(8,i)
1929 DDold(1,2,i)= b(8,i)
1930 DDold(1,1,-i)= b(6,i)
1931 DDold(2,2,-i)=-b(6,i)
1932 DDold(2,1,-i)=-b(8,i)
1933 DDold(1,2,-i)=-b(8,i)
1938 !c Dtilde(1,1,i)= DD(1,1,i)
1939 !c Dtilde(1,2,i)= DD(1,2,i)
1940 !c Dtilde(2,1,i)=-DD(2,1,i)
1941 !c Dtilde(2,2,i)=-DD(2,2,i)
1943 !c Dtilde(1,1,i)=0.0d0
1944 !c Dtilde(1,2,i)=0.0d0
1945 !c Dtilde(2,1,i)=0.0d0
1946 !c Dtilde(2,2,i)=0.0d0
1947 EEold(1,1,i)= b(10,i)+b(11,i)
1948 EEold(2,2,i)=-b(10,i)+b(11,i)
1949 EEold(2,1,i)= b(12,i)-b(13,i)
1950 EEold(1,2,i)= b(12,i)+b(13,i)
1951 EEold(1,1,-i)= b(10,i)+b(11,i)
1952 EEold(2,2,-i)=-b(10,i)+b(11,i)
1953 EEold(2,1,-i)=-b(12,i)+b(13,i)
1954 EEold(1,2,-i)=-b(12,i)-b(13,i)
1955 write(iout,*) "TU DOCHODZE"
1961 !c ee(2,1,i)=ee(1,2,i)
1966 "Coefficients of the cumulants (independent of valence angles)"
1967 do i=-nloctyp+1,nloctyp-1
1968 write (iout,*) 'Type ',onelet(iloctyp(i))
1970 write(iout,'(2f10.5)') B(3,i),B(5,i)
1972 write(iout,'(2f10.5)') B(2,i),B(4,i)
1975 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1979 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1983 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1992 ! Read torsional parameters in old format
1994 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1996 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1997 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1998 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2000 !el from energy module--------
2001 allocate(v1(nterm_old,ntortyp,ntortyp))
2002 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2003 !el---------------------------
2008 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2014 write (iout,'(/a/)') 'Torsional constants:'
2017 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2018 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2024 ! Read torsional parameters
2026 IF (TOR_MODE.eq.0) THEN
2027 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2028 read (itorp,*,end=113,err=113) ntortyp
2029 !el from energy module---------
2030 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2031 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2033 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2034 allocate(vlor2(maxlor,ntortyp,ntortyp))
2035 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2036 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2038 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2039 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2040 !el---------------------------
2043 !el---------------------------
2045 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2047 itortyp(i)=-itortyp(-i)
2049 itortyp(ntyp1)=ntortyp
2050 itortyp(-ntyp1)=-ntortyp
2052 write (iout,*) 'ntortyp',ntortyp
2054 do j=-ntortyp+1,ntortyp-1
2055 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2057 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2058 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2061 do k=1,nterm(i,j,iblock)
2062 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2064 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2065 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2066 v0ij=v0ij+si*v1(k,i,j,iblock)
2068 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2069 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2070 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2072 do k=1,nlor(i,j,iblock)
2073 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2074 vlor2(k,i,j),vlor3(k,i,j)
2075 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2078 v0(-i,-j,iblock)=v0ij
2084 write (iout,'(/a/)') 'Torsional constants:'
2086 do i=-ntortyp,ntortyp
2087 do j=-ntortyp,ntortyp
2088 write (iout,*) 'ityp',i,' jtyp',j
2089 write (iout,*) 'Fourier constants'
2090 do k=1,nterm(i,j,iblock)
2091 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2094 write (iout,*) 'Lorenz constants'
2095 do k=1,nlor(i,j,iblock)
2096 write (iout,'(3(1pe15.5))') &
2097 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2103 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2105 ! 6/23/01 Read parameters for double torsionals
2107 !el from energy module------------
2108 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2109 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2110 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2112 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2113 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2114 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2115 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2116 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2117 !---------------------------------
2121 do j=-ntortyp+1,ntortyp-1
2122 do k=-ntortyp+1,ntortyp-1
2123 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2124 ! write (iout,*) "OK onelett",
2127 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2128 .or. t3.ne.toronelet(k)) then
2129 write (iout,*) "Error in double torsional parameter file",&
2132 call MPI_Finalize(Ierror)
2134 stop "Error in double torsional parameter file"
2136 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2137 ntermd_2(i,j,k,iblock)
2138 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2139 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2140 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2141 ntermd_1(i,j,k,iblock))
2142 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2143 ntermd_1(i,j,k,iblock))
2144 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2145 ntermd_1(i,j,k,iblock))
2146 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2147 ntermd_1(i,j,k,iblock))
2148 ! Martix of D parameters for one dimesional foureir series
2149 do l=1,ntermd_1(i,j,k,iblock)
2150 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2151 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2152 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2153 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2154 ! write(iout,*) "whcodze" ,
2155 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2157 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2158 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2159 v2s(m,l,i,j,k,iblock),&
2160 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2161 ! Martix of D parameters for two dimesional fourier series
2162 do l=1,ntermd_2(i,j,k,iblock)
2164 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2165 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2166 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2167 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2176 write (iout,*) 'Constants for double torsionals'
2179 do j=-ntortyp+1,ntortyp-1
2180 do k=-ntortyp+1,ntortyp-1
2181 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2182 ' nsingle',ntermd_1(i,j,k,iblock),&
2183 ' ndouble',ntermd_2(i,j,k,iblock)
2185 write (iout,*) 'Single angles:'
2186 do l=1,ntermd_1(i,j,k,iblock)
2187 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2188 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2189 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2190 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2193 write (iout,*) 'Pairs of angles:'
2194 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2195 do l=1,ntermd_2(i,j,k,iblock)
2196 write (iout,'(i5,20f10.5)') &
2197 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2200 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2201 do l=1,ntermd_2(i,j,k,iblock)
2202 write (iout,'(i5,20f10.5)') &
2203 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2204 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2214 itype2loc(i)=itortyp(i)
2218 ELSE IF (TOR_MODE.eq.1) THEN
2220 !C read valence-torsional parameters
2221 read (itorp,*,end=121,err=121) ntortyp
2223 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2225 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2226 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2229 itype2loc(i)=itortyp(i)
2233 itortyp(i)=-itortyp(-i)
2235 do i=-ntortyp+1,ntortyp-1
2236 do j=-ntortyp+1,ntortyp-1
2237 !C first we read the cos and sin gamma parameters
2238 read (itorp,'(13x,a)',end=121,err=121) string
2239 write (iout,*) i,j,string
2240 read (itorp,*,end=121,err=121) &
2241 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2242 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2243 do k=1,nterm_kcc(j,i)
2244 do l=1,nterm_kcc_Tb(j,i)
2245 do ll=1,nterm_kcc_Tb(j,i)
2246 read (itorp,*,end=121,err=121) ii,jj,kk, &
2247 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2255 !c AL 4/8/16: Calculate coefficients from one-body parameters
2257 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2258 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2259 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2260 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2261 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2264 print *,i,itortyp(i)
2265 itortyp(i)=itype2loc(i)
2268 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2270 do i=-ntortyp+1,ntortyp-1
2271 do j=-ntortyp+1,ntortyp-1
2274 do k=1,nterm_kcc_Tb(j,i)
2275 do l=1,nterm_kcc_Tb(j,i)
2276 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2277 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2278 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2279 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2282 do k=1,nterm_kcc_Tb(j,i)
2283 do l=1,nterm_kcc_Tb(j,i)
2285 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2286 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2287 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2288 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2290 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2291 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2292 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2293 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2297 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2301 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2306 if (tor_mode.gt.0 .and. lprint) then
2307 !c Print valence-torsional parameters
2308 write (iout,'(a)') &
2309 "Parameters of the valence-torsional potentials"
2310 do i=-ntortyp+1,ntortyp-1
2311 do j=-ntortyp+1,ntortyp-1
2312 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2313 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2314 do k=1,nterm_kcc(j,i)
2315 do l=1,nterm_kcc_Tb(j,i)
2316 do ll=1,nterm_kcc_Tb(j,i)
2317 write (iout,'(3i5,2f15.4)')&
2318 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2326 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2327 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2328 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2329 !el from energy module---------
2330 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2331 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2333 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2334 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2335 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2336 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2338 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2339 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2340 !el---------------------------
2343 !el--------------------
2344 read (itorp_nucl,*,end=113,err=113) &
2345 (itortyp_nucl(i),i=1,ntyp_molec(2))
2346 ! print *,itortyp_nucl(:)
2347 !c write (iout,*) 'ntortyp',ntortyp
2350 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2351 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2354 do k=1,nterm_nucl(i,j)
2355 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2356 v0ij=v0ij+si*v1_nucl(k,i,j)
2359 do k=1,nlor_nucl(i,j)
2360 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2361 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2362 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2368 ! Read of Side-chain backbone correlation parameters
2369 ! Modified 11 May 2012 by Adasko
2372 read (isccor,*,end=119,err=119) nsccortyp
2374 !el from module energy-------------
2375 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2376 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2377 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2378 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2379 !-----------------------------------
2381 !el from module energy-------------
2382 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2384 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2386 isccortyp(i)=-isccortyp(-i)
2388 iscprol=isccortyp(20)
2389 ! write (iout,*) 'ntortyp',ntortyp
2391 !c maxinter is maximum interaction sites
2392 !el from module energy---------
2393 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2394 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2395 -nsccortyp:nsccortyp))
2396 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2397 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2398 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2399 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2400 !-----------------------------------
2404 read (isccor,*,end=119,err=119) &
2405 nterm_sccor(i,j),nlor_sccor(i,j)
2411 nterm_sccor(-i,j)=nterm_sccor(i,j)
2412 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2413 nterm_sccor(i,-j)=nterm_sccor(i,j)
2414 do k=1,nterm_sccor(i,j)
2415 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2417 if (j.eq.iscprol) then
2418 if (i.eq.isccortyp(10)) then
2419 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2420 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2423 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2424 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2425 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2426 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2427 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2428 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2429 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2432 if (i.eq.isccortyp(10)) then
2433 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2434 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2436 if (j.eq.isccortyp(10)) then
2437 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2438 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2440 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2441 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2442 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2443 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2444 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2445 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2449 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2450 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2451 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2452 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2455 do k=1,nlor_sccor(i,j)
2456 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2457 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2458 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2459 (1+vlor3sccor(k,i,j)**2)
2461 v0sccor(l,i,j)=v0ijsccor
2462 v0sccor(l,-i,j)=v0ijsccor1
2463 v0sccor(l,i,-j)=v0ijsccor2
2464 v0sccor(l,-i,-j)=v0ijsccor3
2470 !el from module energy-------------
2471 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2473 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2474 ! write (iout,*) 'ntortyp',ntortyp
2476 !c maxinter is maximum interaction sites
2477 !el from module energy---------
2478 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2479 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2480 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2481 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2482 !-----------------------------------
2486 read (isccor,*,end=119,err=119) &
2487 nterm_sccor(i,j),nlor_sccor(i,j)
2491 do k=1,nterm_sccor(i,j)
2492 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2494 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2497 do k=1,nlor_sccor(i,j)
2498 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2499 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2500 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2501 (1+vlor3sccor(k,i,j)**2)
2503 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2512 write (iout,'(/a/)') 'Torsional constants:'
2515 write (iout,*) 'ityp',i,' jtyp',j
2516 write (iout,*) 'Fourier constants'
2517 do k=1,nterm_sccor(i,j)
2518 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2520 write (iout,*) 'Lorenz constants'
2521 do k=1,nlor_sccor(i,j)
2522 write (iout,'(3(1pe15.5))') &
2523 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2530 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2531 ! interaction energy of the Gly, Ala, and Pro prototypes.
2534 ! Read electrostatic-interaction parameters
2539 write (iout,'(/a)') 'Electrostatic interaction constants:'
2540 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2541 'IT','JT','APP','BPP','AEL6','AEL3'
2543 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2544 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2545 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2546 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2551 app (i,j)=epp(i,j)*rri*rri
2552 bpp (i,j)=-2.0D0*epp(i,j)*rri
2553 ael6(i,j)=elpp6(i,j)*4.2D0**6
2554 ael3(i,j)=elpp3(i,j)*4.2D0**3
2556 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2562 ! Read side-chain interaction parameters.
2564 !el from module energy - COMMON.INTERACT-------
2565 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2566 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2567 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2569 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2570 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2571 allocate(epslip(ntyp,ntyp))
2579 !--------------------------------
2581 read (isidep,*,end=117,err=117) ipot,expon
2582 if (ipot.lt.1 .or. ipot.gt.5) then
2583 write (iout,'(2a)') 'Error while reading SC interaction',&
2584 'potential file - unknown potential type.'
2586 call MPI_Finalize(Ierror)
2592 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2593 ', exponents are ',expon,2*expon
2594 ! goto (10,20,30,30,40) ipot
2596 !----------------------- LJ potential ---------------------------------
2598 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2599 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2600 (sigma0(i),i=1,ntyp)
2602 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2603 write (iout,'(a/)') 'The epsilon array:'
2604 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2605 write (iout,'(/a)') 'One-body parameters:'
2606 write (iout,'(a,4x,a)') 'residue','sigma'
2607 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2610 !----------------------- LJK potential --------------------------------
2612 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2613 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2614 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2616 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2617 write (iout,'(a/)') 'The epsilon array:'
2618 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2619 write (iout,'(/a)') 'One-body parameters:'
2620 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2621 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2625 !---------------------- GB or BP potential -----------------------------
2628 ! print *,"I AM in SCELE",scelemode
2629 if (scelemode.eq.0) then
2631 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2633 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2634 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2635 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2636 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2638 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2641 ! For the GB potential convert sigma'**2 into chi'
2644 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2648 write (iout,'(/a/)') 'Parameters of the BP potential:'
2649 write (iout,'(a/)') 'The epsilon array:'
2650 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2651 write (iout,'(/a)') 'One-body parameters:'
2652 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2654 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2655 chip(i),alp(i),i=1,ntyp)
2658 ! print *,ntyp,"NTYP"
2659 allocate(icharge(ntyp1))
2660 ! print *,ntyp,icharge(i)
2662 read (isidep,*) (icharge(i),i=1,ntyp)
2663 print *,ntyp,icharge(i)
2664 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2665 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2666 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2667 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2668 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2669 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2670 allocate(epsintab(ntyp,ntyp))
2671 allocate(dtail(2,ntyp,ntyp))
2672 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2673 allocate(wqdip(2,ntyp,ntyp))
2674 allocate(wstate(4,ntyp,ntyp))
2675 allocate(dhead(2,2,ntyp,ntyp))
2676 allocate(nstate(ntyp,ntyp))
2677 allocate(debaykap(ntyp,ntyp))
2679 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2680 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2684 ! write (*,*) "Im in ALAB", i, " ", j
2686 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2687 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2688 chis(i,j),chis(j,i), & !2 w tej linii
2689 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2690 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2691 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2692 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2693 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2694 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2695 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2696 IF ((LaTeX).and.(i.gt.24)) then
2697 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2698 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2699 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2700 chis(i,j),chis(j,i) !2 w tej linii
2702 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2707 IF ((LaTeX).and.(i.gt.24)) then
2708 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2709 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2710 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2711 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2712 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2713 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2714 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2721 sigma(i,j) = sigma(j,i)
2722 sigmap1(i,j)=sigmap1(j,i)
2723 sigmap2(i,j)=sigmap2(j,i)
2724 sigiso1(i,j)=sigiso1(j,i)
2725 sigiso2(i,j)=sigiso2(j,i)
2726 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2727 nstate(i,j) = nstate(j,i)
2728 dtail(1,i,j) = dtail(1,j,i)
2729 dtail(2,i,j) = dtail(2,j,i)
2731 alphasur(k,i,j) = alphasur(k,j,i)
2732 wstate(k,i,j) = wstate(k,j,i)
2733 alphiso(k,i,j) = alphiso(k,j,i)
2736 dhead(2,1,i,j) = dhead(1,1,j,i)
2737 dhead(2,2,i,j) = dhead(1,2,j,i)
2738 dhead(1,1,i,j) = dhead(2,1,j,i)
2739 dhead(1,2,i,j) = dhead(2,2,j,i)
2741 epshead(i,j) = epshead(j,i)
2742 sig0head(i,j) = sig0head(j,i)
2745 wqdip(k,i,j) = wqdip(k,j,i)
2748 wquad(i,j) = wquad(j,i)
2749 epsintab(i,j) = epsintab(j,i)
2750 debaykap(i,j)=debaykap(j,i)
2751 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2758 !--------------------- GBV potential -----------------------------------
2760 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2761 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2762 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2763 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2765 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2766 write (iout,'(a/)') 'The epsilon array:'
2767 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2768 write (iout,'(/a)') 'One-body parameters:'
2769 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2770 's||/s_|_^2',' chip ',' alph '
2771 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2772 sigii(i),chip(i),alp(i),i=1,ntyp)
2775 write(iout,*)"Wrong ipot"
2781 !-----------------------------------------------------------------------
2782 ! Calculate the "working" parameters of SC interactions.
2784 !el from module energy - COMMON.INTERACT-------
2785 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2786 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2787 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2788 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2789 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2790 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2792 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2798 if (scelemode.eq.0) then
2807 sc_aa_tube_par(:)=0.0d0
2808 sc_bb_tube_par(:)=0.0d0
2810 !--------------------------------
2815 epslip(i,j)=epslip(j,i)
2818 if (scelemode.eq.0) then
2821 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2822 sigma(j,i)=sigma(i,j)
2823 rs0(i,j)=dwa16*sigma(i,j)
2828 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2829 'Working parameters of the SC interactions:',&
2830 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2835 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2837 ! print *,"SIGMA ZLA?",sigma(i,j)
2845 sigeps=dsign(1.0D0,epsij)
2847 aa_aq(i,j)=epsij*rrij*rrij
2848 ! print *,"ADASKO",epsij,rrij,expon
2849 bb_aq(i,j)=-sigeps*epsij*rrij
2850 aa_aq(j,i)=aa_aq(i,j)
2851 bb_aq(j,i)=bb_aq(i,j)
2852 epsijlip=epslip(i,j)
2853 sigeps=dsign(1.0D0,epsijlip)
2854 epsijlip=dabs(epsijlip)
2855 aa_lip(i,j)=epsijlip*rrij*rrij
2856 bb_lip(i,j)=-sigeps*epsijlip*rrij
2857 aa_lip(j,i)=aa_lip(i,j)
2858 bb_lip(j,i)=bb_lip(i,j)
2859 !C write(iout,*) aa_lip
2860 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2861 sigt1sq=sigma0(i)**2
2862 sigt2sq=sigma0(j)**2
2865 ratsig1=sigt2sq/sigt1sq
2866 ratsig2=1.0D0/ratsig1
2867 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2868 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2869 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2873 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2874 sigmaii(i,j)=rsum_max
2875 sigmaii(j,i)=rsum_max
2877 ! sigmaii(i,j)=r0(i,j)
2878 ! sigmaii(j,i)=r0(i,j)
2880 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2881 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2882 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2883 augm(i,j)=epsij*r_augm**(2*expon)
2884 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2891 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2892 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2893 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2898 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2899 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2900 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2901 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2902 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2903 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2904 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2905 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2906 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2907 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2908 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2909 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2910 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2911 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2912 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2913 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2922 read (isidep_nucl,*) ipot_nucl
2923 ! print *,"TU?!",ipot_nucl
2924 if (ipot_nucl.eq.1) then
2925 do i=1,ntyp_molec(2)
2926 do j=i,ntyp_molec(2)
2927 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2928 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2932 do i=1,ntyp_molec(2)
2933 do j=i,ntyp_molec(2)
2934 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2935 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2936 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2940 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2941 do i=1,ntyp_molec(2)
2942 do j=i,ntyp_molec(2)
2943 rrij=sigma_nucl(i,j)
2947 epsij=4*eps_nucl(i,j)
2948 sigeps=dsign(1.0D0,epsij)
2950 aa_nucl(i,j)=epsij*rrij*rrij
2951 bb_nucl(i,j)=-sigeps*epsij*rrij
2952 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2953 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2954 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2955 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2956 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2957 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2960 aa_nucl(i,j)=aa_nucl(j,i)
2961 bb_nucl(i,j)=bb_nucl(j,i)
2962 ael3_nucl(i,j)=ael3_nucl(j,i)
2963 ael6_nucl(i,j)=ael6_nucl(j,i)
2964 ael63_nucl(i,j)=ael63_nucl(j,i)
2965 ael32_nucl(i,j)=ael32_nucl(j,i)
2966 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2967 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2968 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2969 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2970 eps_nucl(i,j)=eps_nucl(j,i)
2971 sigma_nucl(i,j)=sigma_nucl(j,i)
2972 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2976 write(iout,*) "tube param"
2977 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2978 ccavtubpep,dcavtubpep,tubetranenepep
2979 sigmapeptube=sigmapeptube**6
2980 sigeps=dsign(1.0D0,epspeptube)
2981 epspeptube=dabs(epspeptube)
2982 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2983 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2984 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2986 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2987 ccavtub(i),dcavtub(i),tubetranene(i)
2988 sigmasctube=sigmasctube**6
2989 sigeps=dsign(1.0D0,epssctube)
2990 epssctube=dabs(epssctube)
2991 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2992 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2993 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2995 !-----------------READING SC BASE POTENTIALS-----------------------------
2996 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2997 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2998 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2999 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3000 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3001 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3002 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3003 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3004 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3005 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3006 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3007 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3008 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3009 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3010 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3011 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3014 do i=1,ntyp_molec(1)
3015 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3016 write (*,*) "Im in ", i, " ", j
3017 read(isidep_scbase,*) &
3018 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3019 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3020 write(*,*) "eps",eps_scbase(i,j)
3021 read(isidep_scbase,*) &
3022 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3023 chis_scbase(i,j,1),chis_scbase(i,j,2)
3024 read(isidep_scbase,*) &
3025 dhead_scbasei(i,j), &
3026 dhead_scbasej(i,j), &
3027 rborn_scbasei(i,j),rborn_scbasej(i,j)
3028 read(isidep_scbase,*) &
3029 (wdipdip_scbase(k,i,j),k=1,3), &
3030 (wqdip_scbase(k,i,j),k=1,2)
3031 read(isidep_scbase,*) &
3032 alphapol_scbase(i,j), &
3033 epsintab_scbase(i,j)
3036 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3037 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3039 do i=1,ntyp_molec(1)
3040 do j=1,ntyp_molec(2)-1
3041 epsij=eps_scbase(i,j)
3042 rrij=sigma_scbase(i,j)
3047 sigeps=dsign(1.0D0,epsij)
3049 aa_scbase(i,j)=epsij*rrij*rrij
3050 bb_scbase(i,j)=-sigeps*epsij*rrij
3053 !-----------------READING PEP BASE POTENTIALS-------------------
3054 allocate(eps_pepbase(ntyp_molec(2)))
3055 allocate(sigma_pepbase(ntyp_molec(2)))
3056 allocate(chi_pepbase(ntyp_molec(2),2))
3057 allocate(chipp_pepbase(ntyp_molec(2),2))
3058 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3059 allocate(sigmap1_pepbase(ntyp_molec(2)))
3060 allocate(sigmap2_pepbase(ntyp_molec(2)))
3061 allocate(chis_pepbase(ntyp_molec(2),2))
3062 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3065 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3066 write (*,*) "Im in ", i, " ", j
3067 read(isidep_pepbase,*) &
3068 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3069 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3070 write(*,*) "eps",eps_pepbase(j)
3071 read(isidep_pepbase,*) &
3072 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3073 chis_pepbase(j,1),chis_pepbase(j,2)
3074 read(isidep_pepbase,*) &
3075 (wdipdip_pepbase(k,j),k=1,3)
3077 allocate(aa_pepbase(ntyp_molec(2)))
3078 allocate(bb_pepbase(ntyp_molec(2)))
3080 do j=1,ntyp_molec(2)-1
3081 epsij=eps_pepbase(j)
3082 rrij=sigma_pepbase(j)
3087 sigeps=dsign(1.0D0,epsij)
3089 aa_pepbase(j)=epsij*rrij*rrij
3090 bb_pepbase(j)=-sigeps*epsij*rrij
3092 !--------------READING SC PHOSPHATE-------------------------------------
3093 allocate(eps_scpho(ntyp_molec(1)))
3094 allocate(sigma_scpho(ntyp_molec(1)))
3095 allocate(chi_scpho(ntyp_molec(1),2))
3096 allocate(chipp_scpho(ntyp_molec(1),2))
3097 allocate(alphasur_scpho(4,ntyp_molec(1)))
3098 allocate(sigmap1_scpho(ntyp_molec(1)))
3099 allocate(sigmap2_scpho(ntyp_molec(1)))
3100 allocate(chis_scpho(ntyp_molec(1),2))
3101 allocate(wqq_scpho(ntyp_molec(1)))
3102 allocate(wqdip_scpho(2,ntyp_molec(1)))
3103 allocate(alphapol_scpho(ntyp_molec(1)))
3104 allocate(epsintab_scpho(ntyp_molec(1)))
3105 allocate(dhead_scphoi(ntyp_molec(1)))
3106 allocate(rborn_scphoi(ntyp_molec(1)))
3107 allocate(rborn_scphoj(ntyp_molec(1)))
3108 allocate(alphi_scpho(ntyp_molec(1)))
3112 do j=1,ntyp_molec(1) ! without U then we will take T for U
3113 write (*,*) "Im in scpho ", i, " ", j
3114 read(isidep_scpho,*) &
3115 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3116 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3117 write(*,*) "eps",eps_scpho(j)
3118 read(isidep_scpho,*) &
3119 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3120 chis_scpho(j,1),chis_scpho(j,2)
3121 read(isidep_scpho,*) &
3122 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3123 read(isidep_scpho,*) &
3124 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3128 allocate(aa_scpho(ntyp_molec(1)))
3129 allocate(bb_scpho(ntyp_molec(1)))
3131 do j=1,ntyp_molec(1)
3138 sigeps=dsign(1.0D0,epsij)
3140 aa_scpho(j)=epsij*rrij*rrij
3141 bb_scpho(j)=-sigeps*epsij*rrij
3145 read(isidep_peppho,*) &
3146 eps_peppho,sigma_peppho
3147 read(isidep_peppho,*) &
3148 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3149 read(isidep_peppho,*) &
3150 (wqdip_peppho(k),k=1,2)
3158 sigeps=dsign(1.0D0,epsij)
3160 aa_peppho=epsij*rrij*rrij
3161 bb_peppho=-sigeps*epsij*rrij
3164 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3169 ! Define the SC-p interaction constants (hard-coded; old style)
3172 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3174 ! aad(i,1)=0.3D0*4.0D0**12
3175 ! Following line for constants currently implemented
3176 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3177 aad(i,1)=1.5D0*4.0D0**12
3178 ! aad(i,1)=0.17D0*5.6D0**12
3180 ! "Soft" SC-p repulsion
3182 ! Following line for constants currently implemented
3183 ! aad(i,1)=0.3D0*4.0D0**6
3184 ! "Hard" SC-p repulsion
3185 bad(i,1)=3.0D0*4.0D0**6
3186 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3195 ! 8/9/01 Read the SC-p interaction constants from file
3198 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3201 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3202 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3203 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3204 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3208 write (iout,*) "Parameters of SC-p interactions:"
3210 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3211 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3216 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3218 do i=1,ntyp_molec(2)
3219 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3221 do i=1,ntyp_molec(2)
3222 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3223 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3225 r0pp=1.12246204830937298142*5.16158
3231 ! Define the constants of the disulfide bridge
3235 ! Old arbitrary potential - commented out.
3240 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3241 ! energy surface of diethyl disulfide.
3242 ! A. Liwo and U. Kozlowska, 11/24/03
3260 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3261 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3262 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3263 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3264 allocate(epsintabcat(ntyp,ntyp))
3265 allocate(dtailcat(2,ntyp,ntyp))
3266 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3267 allocate(wqdipcat(2,ntyp,ntyp))
3268 allocate(wstatecat(4,ntyp,ntyp))
3269 allocate(dheadcat(2,2,ntyp,ntyp))
3270 allocate(nstatecat(ntyp,ntyp))
3271 allocate(debaykapcat(ntyp,ntyp))
3273 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3274 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3275 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3276 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3277 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3280 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3281 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3282 if (oldion.eq.0) then
3283 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3284 allocate(icharge(1:ntyp1))
3285 read(iion,*) (icharge(i),i=1,ntyp)
3290 do i=1,ntyp_molec(5)
3291 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3292 print *,msc(i,5),restok(i,5)
3296 do j=1,ntyp_molec(5)
3298 ! do j=1,ntyp_molec(5)
3299 ! write (*,*) "Im in ALAB", i, " ", j
3301 epscat(i,j),sigmacat(i,j), &
3302 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3303 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3305 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3306 ! chiscat(i,j),chiscat(j,i), &
3307 chis1cat(i,j),chis2cat(i,j), &
3309 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3310 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3311 dtailcat(1,i,j),dtailcat(2,i,j), &
3312 epsheadcat(i,j),sig0headcat(i,j), &
3314 ! rborncat(i,j),rborncat(j,i),&
3315 rborn1cat(i,j),rborn2cat(i,j),&
3316 (wqdipcat(k,i,j),k=1,2), &
3317 alphapolcat(i,j),alphapolcat(j,i), &
3318 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3319 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3321 ! write (iout,*) 'i= ', i, ' j= ', j
3322 ! write (iout,*) 'epsi0= ', epscat(1,j)
3323 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3324 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3325 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3326 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3327 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3328 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3329 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3330 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3331 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3332 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3333 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3334 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3335 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3336 ! write (iout,*) 'a1= ', rborncat(j,1)
3337 ! write (iout,*) 'a2= ', rborncat(1,j)
3338 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3339 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3340 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3341 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3345 ! If ((i.eq.1).and.(j.eq.27)) then
3346 ! write (iout,*) 'SEP'
3347 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3348 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3353 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3355 do j=1,ntyp_molec(5)
3359 sigeps=dsign(1.0D0,epsij)
3361 aa_aq_cat(i,j)=epsij*rrij*rrij
3362 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3367 do j=1,ntyp_molec(5)
3369 write (iout,*) 'i= ', i, ' j= ', j
3370 write (iout,*) 'epsi0= ', epscat(i,j)
3371 write (iout,*) 'sigma0= ', sigmacat(i,j)
3372 write (iout,*) 'chi1= ', chi1cat(i,j)
3373 write (iout,*) 'chi1= ', chi2cat(i,j)
3374 write (iout,*) 'chip1= ', chipp1cat(1,j)
3375 write (iout,*) 'chip2= ', chipp2cat(1,j)
3376 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3377 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3378 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3379 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3380 write (iout,*) 'sig1= ', sigmap1cat(1,j)
3381 write (iout,*) 'sig2= ', sigmap2cat(1,j)
3382 write (iout,*) 'chis1= ', chis1cat(1,j)
3383 write (iout,*) 'chis1= ', chis2cat(1,j)
3384 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
3385 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
3386 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3387 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
3388 write (iout,*) 'a1= ', rborn1cat(i,j)
3389 write (iout,*) 'a2= ', rborn2cat(i,j)
3390 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3391 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3392 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
3393 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3394 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3395 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
3398 If ((i.eq.1).and.(j.eq.27)) then
3399 write (iout,*) 'SEP'
3400 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3401 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3411 write (iout,'(/a)') "Disulfide bridge parameters:"
3412 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3413 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3414 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3415 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3418 if (shield_mode.gt.0) then
3419 pi=4.0D0*datan(1.0D0)
3420 !C VSolvSphere the volume of solving sphere
3422 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3423 !C there will be no distinction between proline peptide group and normal peptide
3424 !C group in case of shielding parameters
3425 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3426 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3427 write (iout,*) VSolvSphere,VSolvSphere_div
3428 !C long axis of side chain
3430 long_r_sidechain(i)=vbldsc0(1,i)
3431 ! if (scelemode.eq.0) then
3432 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3433 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3435 ! short_r_sidechain(i)=sigma(i,i)
3437 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3444 111 write (iout,*) "Error reading bending energy parameters."
3446 112 write (iout,*) "Error reading rotamer energy parameters."
3448 113 write (iout,*) "Error reading torsional energy parameters."
3450 114 write (iout,*) "Error reading double torsional energy parameters."
3452 115 write (iout,*) &
3453 "Error reading cumulant (multibody energy) parameters."
3455 116 write (iout,*) "Error reading electrostatic energy parameters."
3457 117 write (iout,*) "Error reading side chain interaction parameters."
3459 118 write (iout,*) "Error reading SCp interaction parameters."
3461 119 write (iout,*) "Error reading SCCOR parameters"
3463 121 write (iout,*) "Error in Czybyshev parameters"
3466 call MPI_Finalize(Ierror)
3470 end subroutine parmread
3472 !-----------------------------------------------------------------------------
3474 !-----------------------------------------------------------------------------
3475 subroutine printmat(ldim,m,n,iout,key,a)
3478 character(len=3),dimension(n) :: key
3479 real(kind=8),dimension(ldim,n) :: a
3481 integer :: i,j,k,m,iout,nlim
3485 write (iout,1000) (key(k),k=i,nlim)
3487 1000 format (/5x,8(6x,a3))
3488 1020 format (/80(1h-)/)
3490 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3493 1010 format (a3,2x,8(f9.4))
3495 end subroutine printmat
3496 !-----------------------------------------------------------------------------
3498 !-----------------------------------------------------------------------------
3500 ! Read the PDB file and convert the peptide geometry into virtual-chain
3503 use energy_data, only: itype,istype
3507 ! use control, only: rescode,sugarcode
3508 ! implicit real*8 (a-h,o-z)
3509 ! include 'DIMENSIONS'
3510 ! include 'COMMON.LOCAL'
3511 ! include 'COMMON.VAR'
3512 ! include 'COMMON.CHAIN'
3513 ! include 'COMMON.INTERACT'
3514 ! include 'COMMON.IOUNITS'
3515 ! include 'COMMON.GEO'
3516 ! include 'COMMON.NAMES'
3517 ! include 'COMMON.CONTROL'
3518 ! include 'COMMON.DISTFIT'
3519 ! include 'COMMON.SETUP'
3520 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3522 logical :: lprn=.true.,fail
3523 real(kind=8),dimension(3) :: e1,e2,e3
3524 real(kind=8) :: dcj,efree_temp
3525 character(len=3) :: seq,res,res2
3526 character(len=5) :: atom
3527 character(len=80) :: card
3528 real(kind=8),dimension(3,20) :: sccor
3529 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3530 integer :: isugar,molecprev,firstion
3531 character*1 :: sugar
3533 real(kind=8),dimension(3) :: ccc
3535 integer,dimension(2,maxres/3) :: hfrag_alloc
3536 integer,dimension(4,maxres/3) :: bfrag_alloc
3537 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3538 real(kind=8),dimension(:,:), allocatable :: c_temporary
3539 integer,dimension(:,:) , allocatable :: itype_temporary
3540 integer,dimension(:),allocatable :: istype_temp
3547 ! write (2,*) "UNRES_PDB",unres_pdb
3567 !-----------------------------
3568 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3569 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3570 if(.not. allocated(istype)) allocate(istype(maxres))
3572 read (ipdbin,'(a80)',end=10) card
3573 write (iout,'(a)') card
3574 if (card(:5).eq.'HELIX') then
3577 read(card(22:25),*) hfrag(1,nhfrag)
3578 read(card(34:37),*) hfrag(2,nhfrag)
3580 if (card(:5).eq.'SHEET') then
3583 read(card(24:26),*) bfrag(1,nbfrag)
3584 read(card(35:37),*) bfrag(2,nbfrag)
3585 !rc----------------------------------------
3586 !rc to be corrected !!!
3587 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3588 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3589 !rc----------------------------------------
3591 if (card(:3).eq.'END') then
3593 else if (card(:3).eq.'TER') then
3598 itype(ires_old,molecule)=ntyp1_molec(molecule)
3599 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3600 nres_molec(molecule)=nres_molec(molecule)+2
3602 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3605 dc(j,ires)=sccor(j,iii)
3608 call sccenter(ires,iii,sccor)
3614 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3615 ! Fish out the ATOM cards.
3616 ! write(iout,*) 'card',card(1:20)
3617 ! print *,"ATU ",card(1:6), CARD(3:6)
3619 if (index(card(1:4),'ATOM').gt.0) then
3620 read (card(12:16),*) atom
3621 ! write (iout,*) "! ",atom," !",ires
3622 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3623 read (card(23:26),*) ires
3624 read (card(18:20),'(a3)') res
3625 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3626 ! & " ires_old",ires_old
3627 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3628 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3629 if (ires-ishift+ishift1.ne.ires_old) then
3630 ! Calculate the CM of the preceding residue.
3631 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3633 ! write (iout,*) "Calculating sidechain center iii",iii
3636 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3639 call sccenter(ires_old,iii,sccor)
3643 ! Start new residue.
3644 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3647 else if (ibeg.eq.1) then
3648 write (iout,*) "BEG ires",ires
3650 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3653 nres_molec(molecule)=nres_molec(molecule)+1
3655 ires=ires-ishift+ishift1
3657 ! write (iout,*) "ishift",ishift," ires",ires,&
3658 ! " ires_old",ires_old
3660 else if (ibeg.eq.2) then
3662 ishift=-ires_old+ires-1 !!!!!
3663 ishift1=ishift1-1 !!!!!
3664 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3665 ires=ires-ishift+ishift1
3666 ! print *,ires,ishift,ishift1
3670 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3671 ires=ires-ishift+ishift1
3674 ! print *,'atom',ires,atom
3675 if (res.eq.'ACE' .or. res.eq.'NHE') then
3678 if (atom.eq.'CA '.or.atom.eq.'N ') then
3680 itype(ires,molecule)=rescode(ires,res,0,molecule)
3682 ! nres_molec(molecule)=nres_molec(molecule)+1
3686 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3687 ! nres_molec(molecule)=nres_molec(molecule)+1
3688 read (card(19:19),'(a1)') sugar
3689 isugar=sugarcode(sugar,ires)
3690 ! if (ibeg.eq.1) then
3694 ! print *,"ires=",ires,istype(ires)
3700 ires=ires-ishift+ishift1
3702 ! write (iout,*) "ires_old",ires_old," ires",ires
3703 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3706 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3707 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3708 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3709 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3710 ! print *,ires,ishift,ishift1
3711 ! write (iout,*) "backbone ",atom
3713 write (iout,'(2i3,2x,a,3f8.3)') &
3714 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3717 nres_molec(molecule)=nres_molec(molecule)+1
3719 sccor(j,iii)=c(j,ires)
3721 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3722 atom.eq."C2'" .or. atom.eq."C3'" &
3723 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3724 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3725 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3726 ! print *,ires,ishift,ishift1
3730 ! sccor(j,iii)=c(j,ires)
3733 c(j,ires)=c(j,ires)+ccc(j)/5.0
3735 print *,counter,molecule
3736 if (counter.eq.5) then
3738 nres_molec(molecule)=nres_molec(molecule)+1
3741 ! sccor(j,iii)=c(j,ires)
3745 ! print *, "ATOM",atom(1:3)
3746 else if (atom.eq."C5'") then
3747 read (card(19:19),'(a1)') sugar
3748 isugar=sugarcode(sugar,ires)
3753 ! print *,ires,istype(ires)
3757 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3758 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3759 nres_molec(molecule)=nres_molec(molecule)+1
3760 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3764 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3766 else if ((atom.eq."C1'").and.unres_pdb) then
3768 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3769 ! write (*,*) card(23:27),ires,itype(ires,1)
3770 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3771 atom.ne.'N' .and. atom.ne.'C' .and. &
3772 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3773 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3774 .and. atom.ne.'P '.and. &
3775 atom(1:1).ne.'H' .and. &
3776 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3778 ! write (iout,*) "sidechain ",atom
3779 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3780 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3781 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3783 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3786 ! print *,"IONS",ions,card(1:6)
3787 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3788 if (firstion.eq.0) then
3792 dc(j,ires)=sccor(j,iii)
3795 call sccenter(ires,iii,sccor)
3798 read (card(12:16),*) atom
3799 ! print *,"HETATOM", atom
3800 read (card(18:20),'(a3)') res
3801 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3802 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3803 .or.(atom(1:2).eq.'K ')) &
3806 if (molecule.ne.5) molecprev=molecule
3808 nres_molec(molecule)=nres_molec(molecule)+1
3809 print *,"HERE",nres_molec(molecule)
3811 itype(ires,molecule)=rescode(ires,res,0,molecule)
3812 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3816 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3817 if (ires.eq.0) return
3818 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3821 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3823 nres_molec(molecule)=nres_molec(molecule)-2
3824 print *,'I have',nres, nres_molec(:)
3826 do k=1,4 ! ions are without dummy
3827 if (nres_molec(k).eq.0) cycle
3829 ! write (iout,*) i,itype(i,1)
3830 ! if (itype(i,1).eq.ntyp1) then
3831 ! write (iout,*) "dummy",i,itype(i,1)
3833 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3834 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3838 if (itype(i,k).eq.ntyp1_molec(k)) then
3839 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3840 if (itype(i+2,k).eq.0) then
3841 ! print *,"masz sieczke"
3843 if (itype(i+2,j).ne.0) then
3845 itype(i+1,j)=ntyp1_molec(j)
3846 nres_molec(k)=nres_molec(k)-1
3847 nres_molec(j)=nres_molec(j)+1
3853 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3854 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3855 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3856 ! if (unres_pdb) then
3857 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3858 ! print *,i,'tu dochodze'
3859 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3867 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3871 dcj=(c(j,i-2)-c(j,i-3))/2.0
3872 if (dcj.eq.0) dcj=1.23591524223
3877 else !itype(i+1,1).eq.ntyp1
3878 ! if (unres_pdb) then
3879 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3880 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3887 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3888 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3892 dcj=(c(j,i+3)-c(j,i+2))/2.0
3893 if (dcj.eq.0) dcj=1.23591524223
3898 endif !itype(i+1,1).eq.ntyp1
3899 endif !itype.eq.ntyp1
3903 ! Calculate the CM of the last side chain.
3907 dc(j,ires)=sccor(j,iii)
3910 call sccenter(ires,iii,sccor)
3916 ! print *,"molecule",molecule
3917 if ((itype(nres,1).ne.10)) then
3919 if (molecule.eq.5) molecule=molecprev
3920 itype(nres,molecule)=ntyp1_molec(molecule)
3921 nres_molec(molecule)=nres_molec(molecule)+1
3922 ! if (unres_pdb) then
3923 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3924 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3931 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3935 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3936 c(j,nres)=c(j,nres-1)+dcj
3937 c(j,2*nres)=c(j,nres)
3941 ! print *,'I have',nres, nres_molec(:)
3943 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3945 if (nres.ne.nres0) then
3946 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3948 stop "Error nres value in WHAM input"
3951 !---------------------------------
3952 !el reallocate tables
3955 ! hfrag_alloc(j,i)=hfrag(j,i)
3958 ! bfrag_alloc(j,i)=bfrag(j,i)
3964 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3965 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3966 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3967 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3971 ! hfrag(j,i)=hfrag_alloc(j,i)
3976 ! bfrag(j,i)=bfrag_alloc(j,i)
3979 !el end reallocate tables
3980 !---------------------------------
3988 c(j,2*nres)=c(j,nres)
3991 if (itype(1,1).eq.ntyp1) then
3995 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3996 call refsys(2,3,4,e1,e2,e3,fail)
4003 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4004 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4008 dcj=(c(j,4)-c(j,3))/2.0
4014 ! First lets assign correct dummy to correct type of chain
4016 if (itype(1,1).eq.ntyp1) then
4017 if (itype(2,1).eq.0) then
4019 if (itype(2,j).ne.0) then
4021 itype(1,j)=ntyp1_molec(j)
4022 nres_molec(1)=nres_molec(1)-1
4023 nres_molec(j)=nres_molec(j)+1
4030 print *,'I have',nres, nres_molec(:)
4032 ! Copy the coordinates to reference coordinates
4038 ! Calculate internal coordinates.
4040 write (iout,'(/a)') &
4041 "Cartesian coordinates of the reference structure"
4042 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4043 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4045 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4046 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4047 (c(j,ires+nres),j=1,3)
4050 ! znamy już nres wiec mozna alokowac tablice
4051 ! Calculate internal coordinates.
4052 if(me.eq.king.or..not.out1file)then
4053 write (iout,'(a)') &
4054 "Backbone and SC coordinates as read from the PDB"
4056 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4057 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4058 (c(j,nres+ires),j=1,3)
4061 ! NOW LETS ROCK! SORTING
4062 allocate(c_temporary(3,2*nres))
4063 allocate(itype_temporary(nres,5))
4064 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4065 if (.not.allocated(istype)) write(iout,*) &
4066 "SOMETHING WRONG WITH ISTYTPE"
4067 allocate(istype_temp(nres))
4068 itype_temporary(:,:)=0
4072 if (itype(i,k).ne.0) then
4074 c_temporary(j,seqalingbegin)=c(j,i)
4075 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4078 itype_temporary(seqalingbegin,k)=itype(i,k)
4079 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4080 istype_temp(seqalingbegin)=istype(i)
4081 molnum(seqalingbegin)=k
4082 seqalingbegin=seqalingbegin+1
4088 c(j,i)=c_temporary(j,i)
4093 itype(i,k)=itype_temporary(i,k)
4094 istype(i)=istype_temp(i)
4097 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4098 ! I have only ions now dummy atoms in the system
4100 itype(1,5)=itype(2,5)
4103 itype(i,5)=itype(i+1,5)
4107 nres_molec(1)=nres_molec(1)-1
4109 ! if (itype(1,1).eq.ntyp1) then
4112 ! if (unres_pdb) then
4113 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4114 ! call refsys(2,3,4,e1,e2,e3,fail)
4121 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4125 ! dcj=(c(j,4)-c(j,3))/2.0
4127 ! c(j,nres+1)=c(j,1)
4133 write (iout,'(/a)') &
4134 "Cartesian coordinates of the reference structure after sorting"
4135 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4136 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4138 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4139 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4140 (c(j,ires+nres),j=1,3)
4144 ! print *,seqalingbegin,nres
4145 if(.not.allocated(vbld)) then
4146 allocate(vbld(2*nres))
4151 if(.not.allocated(vbld_inv)) then
4152 allocate(vbld_inv(2*nres))
4158 if(.not.allocated(theta)) then
4159 allocate(theta(nres+2))
4163 if(.not.allocated(phi)) allocate(phi(nres+2))
4164 if(.not.allocated(alph)) allocate(alph(nres+2))
4165 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4166 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4167 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4168 if(.not.allocated(costtab)) allocate(costtab(nres))
4169 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4170 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4171 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4172 if(.not.allocated(xxref)) allocate(xxref(nres))
4173 if(.not.allocated(yyref)) allocate(yyref(nres))
4174 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4175 if(.not.allocated(dc_norm)) then
4176 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4177 allocate(dc_norm(3,0:2*nres+2))
4181 call int_from_cart(.true.,.false.)
4182 call sc_loc_geom(.false.)
4184 thetaref(i)=theta(i)
4194 dc(j,i)=c(j,i+1)-c(j,i)
4195 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4200 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4201 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4203 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4207 ! Copy the coordinates to reference coordinates
4208 ! Splits to single chain if occurs
4212 ! cref(j,i,cou)=c(j,i)
4216 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4217 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4218 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4219 !-----------------------------
4223 write (iout,*) "symetr", symetr
4226 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4228 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4231 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4236 cref(j,i,cou)=c(j,i)
4237 cref(j,i+nres,cou)=c(j,i+nres)
4239 chain_rep(j,lll,kkk)=c(j,i)
4240 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4244 write (iout,*) chain_length
4245 if (chain_length.eq.0) chain_length=nres
4247 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4248 chain_rep(j,chain_length+nres,symetr) &
4249 =chain_rep(j,chain_length+nres,1)
4252 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4254 ! do kkk=1,chain_length
4255 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4259 ! makes copy of chains
4260 write (iout,*) "symetr", symetr
4265 if (symetr.gt.1) then
4272 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4278 write (iout,*) i,icha
4279 do lll=1,chain_length
4281 if (cou.le.nres) then
4283 kupa=mod(lll,chain_length)
4284 iprzes=(kkk-1)*chain_length+lll
4285 if (kupa.eq.0) kupa=chain_length
4286 write (iout,*) "kupa", kupa
4287 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4288 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4295 !-koniec robienia kopii
4298 write (iout,*) "nowa struktura", nperm
4300 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4302 cref(3,i,kkk),cref(1,nres+i,kkk),&
4303 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4305 100 format (//' alpha-carbon coordinates ',&
4306 ' centroid coordinates'/ &
4307 ' ', 6X,'X',11X,'Y',11X,'Z', &
4308 10X,'X',11X,'Y',11X,'Z')
4309 110 format (a,'(',i5,')',6f12.5)
4315 bfrag(i,j)=bfrag(i,j)-ishift
4321 hfrag(i,j)=hfrag(i,j)-ishift
4327 end subroutine readpdb
4328 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4329 !-----------------------------------------------------------------------------
4331 !-----------------------------------------------------------------------------
4332 subroutine read_control
4346 use random, only: random_init
4347 ! implicit real*8 (a-h,o-z)
4348 ! include 'DIMENSIONS'
4350 use prng, only:prng_restart
4352 logical :: OKRandom!, prng_restart
4355 ! include 'COMMON.IOUNITS'
4356 ! include 'COMMON.TIME1'
4357 ! include 'COMMON.THREAD'
4358 ! include 'COMMON.SBRIDGE'
4359 ! include 'COMMON.CONTROL'
4360 ! include 'COMMON.MCM'
4361 ! include 'COMMON.MAP'
4362 ! include 'COMMON.HEADER'
4363 ! include 'COMMON.CSA'
4364 ! include 'COMMON.CHAIN'
4365 ! include 'COMMON.MUCA'
4366 ! include 'COMMON.MD'
4367 ! include 'COMMON.FFIELD'
4368 ! include 'COMMON.INTERACT'
4369 ! include 'COMMON.SETUP'
4370 !el integer :: KDIAG,ICORFL,IXDR
4371 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4372 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4373 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4374 ! character(len=80) :: ucase
4375 character(len=640) :: controlcard
4377 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4383 read (INP,'(a)') titel
4384 call card_concat(controlcard,.true.)
4385 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4386 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4387 call reada(controlcard,'SEED',seed,0.0D0)
4388 call random_init(seed)
4389 ! Set up the time limit (caution! The time must be input in minutes!)
4390 read_cart=index(controlcard,'READ_CART').gt.0
4391 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4392 call readi(controlcard,'SYM',symetr,1)
4393 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4394 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4395 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4396 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4397 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4398 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4399 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4400 call reada(controlcard,'DRMS',drms,0.1D0)
4401 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4402 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4403 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4404 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4405 write (iout,'(a,f10.1)')'DRMS = ',drms
4406 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4407 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4409 call readi(controlcard,'NZ_START',nz_start,0)
4410 call readi(controlcard,'NZ_END',nz_end,0)
4411 call readi(controlcard,'IZ_SC',iz_sc,0)
4412 timlim=60.0D0*timlim
4413 safety = 60.0d0*safety
4416 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4417 !C SHIELD keyword sets if the shielding effect of side-chains is used
4418 !C 0 denots no shielding is used all peptide are equally despite the
4419 !C solvent accesible area
4420 !C 1 the newly introduced function
4421 !C 2 reseved for further possible developement
4422 call readi(controlcard,'SHIELD',shield_mode,0)
4423 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4424 write(iout,*) "shield_mode",shield_mode
4425 call readi(controlcard,'TORMODE',tor_mode,0)
4426 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4427 write(iout,*) "torsional and valence angle mode",tor_mode
4429 !C Varibles set size of box
4430 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4431 protein=index(controlcard,"PROTEIN").gt.0
4432 ions=index(controlcard,"IONS").gt.0
4433 call readi(controlcard,'OLDION',oldion,1)
4434 nucleic=index(controlcard,"NUCLEIC").gt.0
4435 write (iout,*) "with_theta_constr ",with_theta_constr
4436 AFMlog=(index(controlcard,'AFM'))
4437 selfguide=(index(controlcard,'SELFGUIDE'))
4438 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4439 call readi(controlcard,'GENCONSTR',genconstr,0)
4440 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4441 call reada(controlcard,'BOXY',boxysize,100.0d0)
4442 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4443 call readi(controlcard,'TUBEMOD',tubemode,0)
4444 print *,"SCELE",scelemode
4445 call readi(controlcard,"SCELEMODE",scelemode,0)
4446 print *,"SCELE",scelemode
4448 ! elemode = 0 is orignal UNRES electrostatics
4449 ! elemode = 1 is "Momo" potentials in progress
4450 ! elemode = 2 is in development EVALD
4453 write (iout,*) TUBEmode,"TUBEMODE"
4454 if (TUBEmode.gt.0) then
4455 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4456 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4457 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4458 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4459 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4460 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4461 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4462 buftubebot=bordtubebot+tubebufthick
4463 buftubetop=bordtubetop-tubebufthick
4466 ! CUTOFFF ON ELECTROSTATICS
4467 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4468 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4469 write(iout,*) "R_CUT_ELE=",r_cut_ele
4470 ! Lipidic parameters
4471 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4472 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4473 if (lipthick.gt.0.0d0) then
4474 bordliptop=(boxzsize+lipthick)/2.0
4475 bordlipbot=bordliptop-lipthick
4476 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4477 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4478 buflipbot=bordlipbot+lipbufthick
4479 bufliptop=bordliptop-lipbufthick
4480 if ((lipbufthick*2.0d0).gt.lipthick) &
4481 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4482 endif !lipthick.gt.0
4483 write(iout,*) "bordliptop=",bordliptop
4484 write(iout,*) "bordlipbot=",bordlipbot
4485 write(iout,*) "bufliptop=",bufliptop
4486 write(iout,*) "buflipbot=",buflipbot
4487 write (iout,*) "SHIELD MODE",shield_mode
4489 !C-------------------------
4490 minim=(index(controlcard,'MINIMIZE').gt.0)
4491 dccart=(index(controlcard,'CART').gt.0)
4492 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4493 overlapsc=.not.overlapsc
4494 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4495 searchsc=.not.searchsc
4496 sideadd=(index(controlcard,'SIDEADD').gt.0)
4497 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4498 outpdb=(index(controlcard,'PDBOUT').gt.0)
4499 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4500 pdbref=(index(controlcard,'PDBREF').gt.0)
4501 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4502 indpdb=index(controlcard,'PDBSTART')
4503 extconf=(index(controlcard,'EXTCONF').gt.0)
4504 call readi(controlcard,'IPRINT',iprint,0)
4505 call readi(controlcard,'MAXGEN',maxgen,10000)
4506 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4507 call readi(controlcard,"KDIAG",kdiag,0)
4508 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4509 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4510 write (iout,*) "RESCALE_MODE",rescale_mode
4511 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4512 if (index(controlcard,'REGULAR').gt.0.0D0) then
4513 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4517 if (index(controlcard,'CHECKGRAD').gt.0) then
4519 if (index(controlcard,'CART').gt.0) then
4521 elseif (index(controlcard,'CARINT').gt.0) then
4526 elseif (index(controlcard,'THREAD').gt.0) then
4528 call readi(controlcard,'THREAD',nthread,0)
4529 if (nthread.gt.0) then
4530 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4533 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4534 stop 'Error termination in Read_Control.'
4536 else if (index(controlcard,'MCMA').gt.0) then
4538 else if (index(controlcard,'MCEE').gt.0) then
4540 else if (index(controlcard,'MULTCONF').gt.0) then
4542 else if (index(controlcard,'MAP').gt.0) then
4544 call readi(controlcard,'MAP',nmap,0)
4545 else if (index(controlcard,'CSA').gt.0) then
4547 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4549 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4552 !fcm else if (index(controlcard,'MCMF').gt.0) then
4554 else if (index(controlcard,'SOFTREG').gt.0) then
4556 else if (index(controlcard,'CHECK_BOND').gt.0) then
4558 else if (index(controlcard,'TEST').gt.0) then
4560 else if (index(controlcard,'MD').gt.0) then
4562 else if (index(controlcard,'RE ').gt.0) then
4566 lmuca=index(controlcard,'MUCA').gt.0
4567 call readi(controlcard,'MUCADYN',mucadyn,0)
4568 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4569 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4571 write (iout,*) 'MUCADYN=',mucadyn
4572 write (iout,*) 'MUCASMOOTH=',muca_smooth
4575 iscode=index(controlcard,'ONE_LETTER')
4576 indphi=index(controlcard,'PHI')
4577 indback=index(controlcard,'BACK')
4578 iranconf=index(controlcard,'RAND_CONF')
4579 i2ndstr=index(controlcard,'USE_SEC_PRED')
4580 gradout=index(controlcard,'GRADOUT').gt.0
4581 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4582 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4583 if (me.eq.king .or. .not.out1file ) &
4584 write (iout,*) "DISTCHAINMAX",distchainmax
4586 if(me.eq.king.or..not.out1file) &
4587 write (iout,'(2a)') diagmeth(kdiag),&
4588 ' routine used to diagonalize matrices.'
4589 if (shield_mode.gt.0) then
4590 pi=4.0D0*datan(1.0D0)
4591 !C VSolvSphere the volume of solving sphere
4593 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4594 !C there will be no distinction between proline peptide group and normal peptide
4595 !C group in case of shielding parameters
4596 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4597 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4598 write (iout,*) VSolvSphere,VSolvSphere_div
4599 !C long axis of side chain
4601 ! long_r_sidechain(i)=vbldsc0(1,i)
4602 ! short_r_sidechain(i)=sigma0(i)
4603 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4610 end subroutine read_control
4611 !-----------------------------------------------------------------------------
4612 subroutine read_REMDpar
4614 ! Read REMD settings
4621 use control_data, only:out1file
4622 ! implicit real*8 (a-h,o-z)
4623 ! include 'DIMENSIONS'
4624 ! include 'COMMON.IOUNITS'
4625 ! include 'COMMON.TIME1'
4626 ! include 'COMMON.MD'
4629 !el include 'COMMON.LANGEVIN'
4631 !el include 'COMMON.LANGEVIN.lang0'
4633 ! include 'COMMON.INTERACT'
4634 ! include 'COMMON.NAMES'
4635 ! include 'COMMON.GEO'
4636 ! include 'COMMON.REMD'
4637 ! include 'COMMON.CONTROL'
4638 ! include 'COMMON.SETUP'
4639 ! character(len=80) :: ucase
4640 character(len=320) :: controlcard
4641 character(len=3200) :: controlcard1
4642 integer :: iremd_m_total
4645 ! real(kind=8) :: var,ene
4647 if(me.eq.king.or..not.out1file) &
4648 write (iout,*) "REMD setup"
4650 call card_concat(controlcard,.true.)
4651 call readi(controlcard,"NREP",nrep,3)
4652 call readi(controlcard,"NSTEX",nstex,1000)
4653 call reada(controlcard,"RETMIN",retmin,10.0d0)
4654 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4655 mremdsync=(index(controlcard,'SYNC').gt.0)
4656 call readi(controlcard,"NSYN",i_sync_step,100)
4657 restart1file=(index(controlcard,'REST1FILE').gt.0)
4658 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4659 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4660 if(max_cache_traj_use.gt.max_cache_traj) &
4661 max_cache_traj_use=max_cache_traj
4662 if(me.eq.king.or..not.out1file) then
4663 !d if (traj1file) then
4664 !rc caching is in testing - NTWX is not ignored
4665 !d write (iout,*) "NTWX value is ignored"
4666 !d write (iout,*) " trajectory is stored to one file by master"
4667 !d write (iout,*) " before exchange at NSTEX intervals"
4669 write (iout,*) "NREP= ",nrep
4670 write (iout,*) "NSTEX= ",nstex
4671 write (iout,*) "SYNC= ",mremdsync
4672 write (iout,*) "NSYN= ",i_sync_step
4673 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4676 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4677 if (index(controlcard,'TLIST').gt.0) then
4679 call card_concat(controlcard1,.true.)
4680 read(controlcard1,*) (remd_t(i),i=1,nrep)
4681 if(me.eq.king.or..not.out1file) &
4682 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4685 if (index(controlcard,'MLIST').gt.0) then
4687 call card_concat(controlcard1,.true.)
4688 read(controlcard1,*) (remd_m(i),i=1,nrep)
4689 if(me.eq.king.or..not.out1file) then
4690 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4693 iremd_m_total=iremd_m_total+remd_m(i)
4695 write (iout,*) 'Total number of replicas ',iremd_m_total
4698 if(me.eq.king.or..not.out1file) &
4699 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4701 end subroutine read_REMDpar
4702 !-----------------------------------------------------------------------------
4703 subroutine read_MDpar
4707 use control_data, only: r_cut,rlamb,out1file
4709 use geometry_data, only: pi
4711 ! implicit real*8 (a-h,o-z)
4712 ! include 'DIMENSIONS'
4713 ! include 'COMMON.IOUNITS'
4714 ! include 'COMMON.TIME1'
4715 ! include 'COMMON.MD'
4718 !el include 'COMMON.LANGEVIN'
4720 !el include 'COMMON.LANGEVIN.lang0'
4722 ! include 'COMMON.INTERACT'
4723 ! include 'COMMON.NAMES'
4724 ! include 'COMMON.GEO'
4725 ! include 'COMMON.SETUP'
4726 ! include 'COMMON.CONTROL'
4727 ! include 'COMMON.SPLITELE'
4728 ! character(len=80) :: ucase
4729 character(len=320) :: controlcard
4734 call card_concat(controlcard,.true.)
4735 call readi(controlcard,"NSTEP",n_timestep,1000000)
4736 call readi(controlcard,"NTWE",ntwe,100)
4737 call readi(controlcard,"NTWX",ntwx,1000)
4738 call reada(controlcard,"DT",d_time,1.0d-1)
4739 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4740 call reada(controlcard,"DAMAX",damax,1.0d1)
4741 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4742 call readi(controlcard,"LANG",lang,0)
4743 RESPA = index(controlcard,"RESPA") .gt. 0
4744 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4745 ntime_split0=ntime_split
4746 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4747 ntime_split0=ntime_split
4748 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4749 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4750 rest = index(controlcard,"REST").gt.0
4751 tbf = index(controlcard,"TBF").gt.0
4752 usampl = index(controlcard,"USAMPL").gt.0
4753 mdpdb = index(controlcard,"MDPDB").gt.0
4754 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4755 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4756 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4757 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4758 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4759 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4760 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4761 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4762 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4763 large = index(controlcard,"LARGE").gt.0
4764 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4765 rattle = index(controlcard,"RATTLE").gt.0
4766 preminim=(index(controlcard,'PREMINIM').gt.0)
4767 write (iout,*) "PREMINIM ",preminim
4768 dccart=(index(controlcard,'CART').gt.0)
4769 if (preminim) call read_minim
4770 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4776 if(me.eq.king.or..not.out1file) then
4778 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4780 write (iout,'(a)') "The units are:"
4781 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4782 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4783 " acceleration: angstrom/(48.9 fs)**2"
4784 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4786 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4787 write (iout,'(a60,f10.5,a)') &
4788 "Initial time step of numerical integration:",d_time,&
4790 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4792 write (iout,'(2a,i4,a)') &
4793 "A-MTS algorithm used; initial time step for fast-varying",&
4794 " short-range forces split into",ntime_split," steps."
4795 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4796 r_cut," lambda",rlamb
4798 write (iout,'(2a,f10.5)') &
4799 "Maximum acceleration threshold to reduce the time step",&
4800 "/increase split number:",damax
4801 write (iout,'(2a,f10.5)') &
4802 "Maximum predicted energy drift to reduce the timestep",&
4803 "/increase split number:",edriftmax
4804 write (iout,'(a60,f10.5)') &
4805 "Maximum velocity threshold to reduce velocities:",dvmax
4806 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4807 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4808 if (rattle) write (iout,'(a60)') &
4809 "Rattle algorithm used to constrain the virtual bonds"
4813 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4814 call reada(controlcard,"RWAT",rwat,1.4d0)
4815 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4816 surfarea=index(controlcard,"SURFAREA").gt.0
4817 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4818 if(me.eq.king.or..not.out1file)then
4819 write (iout,'(/a,$)') "Langevin dynamics calculation"
4821 write (iout,'(a/)') &
4822 " with direct integration of Langevin equations"
4823 else if (lang.eq.2) then
4824 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4825 else if (lang.eq.3) then
4826 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4827 else if (lang.eq.4) then
4828 write (iout,'(a/)') " in overdamped mode"
4830 write (iout,'(//a,i5)') &
4831 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4834 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4835 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4836 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4837 write (iout,'(a60,f10.5)') &
4838 "Scaling factor of the friction forces:",scal_fric
4839 if (surfarea) write (iout,'(2a,i10,a)') &
4840 "Friction coefficients will be scaled by solvent-accessible",&
4841 " surface area every",reset_fricmat," steps."
4843 ! Calculate friction coefficients and bounds of stochastic forces
4844 eta=6*pi*cPoise*etawat
4845 if(me.eq.king.or..not.out1file) &
4846 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4849 do j=1,5 !types of molecules
4850 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4851 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4853 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4854 do j=1,5 !types of molecules
4856 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4857 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4861 if(me.eq.king.or..not.out1file)then
4862 write (iout,'(/2a/)') &
4863 "Radii of site types and friction coefficients and std's of",&
4864 " stochastic forces of fully exposed sites"
4865 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4867 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4868 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4872 if(me.eq.king.or..not.out1file)then
4873 write (iout,'(a)') "Berendsen bath calculation"
4874 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4875 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4877 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4878 count_reset_moment," steps"
4880 write (iout,'(a,i10,a)') &
4881 "Velocities will be reset at random every",count_reset_vel,&
4885 if(me.eq.king.or..not.out1file) &
4886 write (iout,'(a31)') "Microcanonical mode calculation"
4888 if(me.eq.king.or..not.out1file)then
4889 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4891 write(iout,*) "MD running with constraints."
4892 write(iout,*) "Equilibration time ", eq_time, " mtus."
4893 write(iout,*) "Constraining ", nfrag," fragments."
4894 write(iout,*) "Length of each fragment, weight and q0:"
4896 write (iout,*) "Set of restraints #",iset
4898 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4899 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4901 write(iout,*) "constraints between ", npair, "fragments."
4902 write(iout,*) "constraint pairs, weights and q0:"
4904 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4905 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4907 write(iout,*) "angle constraints within ", nfrag_back,&
4908 "backbone fragments."
4909 write(iout,*) "fragment, weights:"
4911 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4912 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4913 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4916 iset=mod(kolor,nset)+1
4919 if(me.eq.king.or..not.out1file) &
4920 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4922 end subroutine read_MDpar
4923 !-----------------------------------------------------------------------------
4927 ! implicit real*8 (a-h,o-z)
4928 ! include 'DIMENSIONS'
4929 ! include 'COMMON.MAP'
4930 ! include 'COMMON.IOUNITS'
4931 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4932 character(len=80) :: mapcard !,ucase
4935 ! real(kind=8) :: var,ene
4938 read (inp,'(a)') mapcard
4939 mapcard=ucase(mapcard)
4940 if (index(mapcard,'PHI').gt.0) then
4942 else if (index(mapcard,'THE').gt.0) then
4944 else if (index(mapcard,'ALP').gt.0) then
4946 else if (index(mapcard,'OME').gt.0) then
4949 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4950 stop 'Error - illegal variable spec in MAP card.'
4952 call readi (mapcard,'RES1',res1(imap),0)
4953 call readi (mapcard,'RES2',res2(imap),0)
4954 if (res1(imap).eq.0) then
4955 res1(imap)=res2(imap)
4956 else if (res2(imap).eq.0) then
4957 res2(imap)=res1(imap)
4959 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4960 write (iout,'(a)') &
4961 'Error - illegal definition of variable group in MAP.'
4962 stop 'Error - illegal definition of variable group in MAP.'
4964 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4965 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4966 call readi(mapcard,'NSTEP',nstep(imap),0)
4967 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4968 write (iout,'(a)') &
4969 'Illegal boundary and/or step size specification in MAP.'
4970 stop 'Illegal boundary and/or step size specification in MAP.'
4974 end subroutine map_read
4975 !-----------------------------------------------------------------------------
4978 use control_data, only: vdisulf
4980 ! implicit real*8 (a-h,o-z)
4981 ! include 'DIMENSIONS'
4982 ! include 'COMMON.IOUNITS'
4983 ! include 'COMMON.GEO'
4984 ! include 'COMMON.CSA'
4985 ! include 'COMMON.BANK'
4986 ! include 'COMMON.CONTROL'
4987 ! character(len=80) :: ucase
4988 character(len=620) :: mcmcard
4990 ! integer :: ntf,ik,iw_pdb
4991 ! real(kind=8) :: var,ene
4993 call card_concat(mcmcard,.true.)
4995 call readi(mcmcard,'NCONF',nconf,50)
4996 call readi(mcmcard,'NADD',nadd,0)
4997 call readi(mcmcard,'JSTART',jstart,1)
4998 call readi(mcmcard,'JEND',jend,1)
4999 call readi(mcmcard,'NSTMAX',nstmax,500000)
5000 call readi(mcmcard,'N0',n0,1)
5001 call readi(mcmcard,'N1',n1,6)
5002 call readi(mcmcard,'N2',n2,4)
5003 call readi(mcmcard,'N3',n3,0)
5004 call readi(mcmcard,'N4',n4,0)
5005 call readi(mcmcard,'N5',n5,0)
5006 call readi(mcmcard,'N6',n6,10)
5007 call readi(mcmcard,'N7',n7,0)
5008 call readi(mcmcard,'N8',n8,0)
5009 call readi(mcmcard,'N9',n9,0)
5010 call readi(mcmcard,'N14',n14,0)
5011 call readi(mcmcard,'N15',n15,0)
5012 call readi(mcmcard,'N16',n16,0)
5013 call readi(mcmcard,'N17',n17,0)
5014 call readi(mcmcard,'N18',n18,0)
5016 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5018 call readi(mcmcard,'NDIFF',ndiff,2)
5019 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5020 call readi(mcmcard,'IS1',is1,1)
5021 call readi(mcmcard,'IS2',is2,8)
5022 call readi(mcmcard,'NRAN0',nran0,4)
5023 call readi(mcmcard,'NRAN1',nran1,2)
5024 call readi(mcmcard,'IRR',irr,1)
5025 call readi(mcmcard,'NSEED',nseed,20)
5026 call readi(mcmcard,'NTOTAL',ntotal,10000)
5027 call reada(mcmcard,'CUT1',cut1,2.0d0)
5028 call reada(mcmcard,'CUT2',cut2,5.0d0)
5029 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5030 call readi(mcmcard,'ICMAX',icmax,3)
5031 call readi(mcmcard,'IRESTART',irestart,0)
5032 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5035 call reada(mcmcard,'DELE',dele,20.0d0)
5036 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5037 call readi(mcmcard,'IREF',iref,0)
5038 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5039 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5040 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5041 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5042 write (iout,*) "NCONF_IN",nconf_in
5044 end subroutine csaread
5045 !-----------------------------------------------------------------------------
5049 use control_data, only: MaxMoveType
5052 ! implicit real*8 (a-h,o-z)
5053 ! include 'DIMENSIONS'
5054 ! include 'COMMON.MCM'
5055 ! include 'COMMON.MCE'
5056 ! include 'COMMON.IOUNITS'
5057 ! character(len=80) :: ucase
5058 character(len=320) :: mcmcard
5061 ! real(kind=8) :: var,ene
5063 call card_concat(mcmcard,.true.)
5064 call readi(mcmcard,'MAXACC',maxacc,100)
5065 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5066 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5067 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5068 call readi(mcmcard,'MAXREPM',maxrepm,200)
5069 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5070 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5071 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5072 call reada(mcmcard,'E_UP',e_up,5.0D0)
5073 call reada(mcmcard,'DELTE',delte,0.1D0)
5074 call readi(mcmcard,'NSWEEP',nsweep,5)
5075 call readi(mcmcard,'NSTEPH',nsteph,0)
5076 call readi(mcmcard,'NSTEPC',nstepc,0)
5077 call reada(mcmcard,'TMIN',tmin,298.0D0)
5078 call reada(mcmcard,'TMAX',tmax,298.0D0)
5079 call readi(mcmcard,'NWINDOW',nwindow,0)
5080 call readi(mcmcard,'PRINT_MC',print_mc,0)
5081 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5082 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5083 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5084 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5085 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5086 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5087 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5088 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5089 if (nwindow.gt.0) then
5090 allocate(winstart(nwindow)) !!el (maxres)
5091 allocate(winend(nwindow)) !!el
5092 allocate(winlen(nwindow)) !!el
5093 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5095 winlen(i)=winend(i)-winstart(i)+1
5098 if (tmax.lt.tmin) tmax=tmin
5099 if (tmax.eq.tmin) then
5103 if (nstepc.gt.0 .and. nsteph.gt.0) then
5104 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5105 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5107 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5108 ! Probabilities of different move types
5109 sumpro_type(0)=0.0D0
5110 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5111 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5112 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5113 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5114 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5115 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5116 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5118 print *,'i',i,' sumprotype',sumpro_type(i)
5119 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5120 print *,'i',i,' sumprotype',sumpro_type(i)
5123 end subroutine mcmread
5124 !-----------------------------------------------------------------------------
5125 subroutine read_minim
5128 ! implicit real*8 (a-h,o-z)
5129 ! include 'DIMENSIONS'
5130 ! include 'COMMON.MINIM'
5131 ! include 'COMMON.IOUNITS'
5132 ! character(len=80) :: ucase
5133 character(len=320) :: minimcard
5135 ! integer :: ntf,ik,iw_pdb
5136 ! real(kind=8) :: var,ene
5138 call card_concat(minimcard,.true.)
5139 call readi(minimcard,'MAXMIN',maxmin,2000)
5140 call readi(minimcard,'MAXFUN',maxfun,5000)
5141 call readi(minimcard,'MINMIN',minmin,maxmin)
5142 call readi(minimcard,'MINFUN',minfun,maxmin)
5143 call reada(minimcard,'TOLF',tolf,1.0D-2)
5144 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5145 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5146 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5147 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5148 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5149 'Options in energy minimization:'
5150 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5151 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5152 'MinMin:',MinMin,' MinFun:',MinFun,&
5153 ' TolF:',TolF,' RTolF:',RTolF
5155 end subroutine read_minim
5156 !-----------------------------------------------------------------------------
5157 subroutine openunits
5159 use MD_data, only: usampl
5162 use control_data, only:out1file
5163 use control, only: getenv_loc
5164 ! implicit real*8 (a-h,o-z)
5165 ! include 'DIMENSIONS'
5168 character(len=16) :: form,nodename
5169 integer :: nodelen,ierror,npos
5171 ! include 'COMMON.SETUP'
5172 ! include 'COMMON.IOUNITS'
5173 ! include 'COMMON.MD'
5174 ! include 'COMMON.CONTROL'
5175 integer :: lenpre,lenpot,lentmp !,ilen
5177 character(len=3) :: out1file_text !,ucase
5178 character(len=3) :: ll
5181 ! integer :: ntf,ik,iw_pdb
5182 ! real(kind=8) :: var,ene
5184 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5185 call getenv_loc("PREFIX",prefix)
5187 call getenv_loc("POT",pot)
5188 call getenv_loc("DIRTMP",tmpdir)
5189 call getenv_loc("CURDIR",curdir)
5190 call getenv_loc("OUT1FILE",out1file_text)
5191 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5192 out1file_text=ucase(out1file_text)
5193 if (out1file_text(1:1).eq."Y") then
5196 out1file=fg_rank.gt.0
5201 if (lentmp.gt.0) then
5202 write (*,'(80(1h!))')
5203 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5204 write (*,'(80(1h!))')
5205 write (*,*)"All output files will be on node /tmp directory."
5207 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5208 if (me.eq.king) then
5209 write (*,*) "The master node is ",nodename
5210 else if (fg_rank.eq.0) then
5211 write (*,*) "I am the CG slave node ",nodename
5213 write (*,*) "I am the FG slave node ",nodename
5216 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5217 lenpre = lentmp+lenpre+1
5219 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5220 ! Get the names and open the input files
5221 #if defined(WINIFL) || defined(WINPGI)
5222 open(1,file=pref_orig(:ilen(pref_orig))// &
5223 '.inp',status='old',readonly,shared)
5224 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5225 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5226 ! Get parameter filenames and open the parameter files.
5227 call getenv_loc('BONDPAR',bondname)
5228 open (ibond,file=bondname,status='old',readonly,shared)
5229 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5230 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5231 call getenv_loc('THETPAR',thetname)
5232 open (ithep,file=thetname,status='old',readonly,shared)
5233 call getenv_loc('ROTPAR',rotname)
5234 open (irotam,file=rotname,status='old',readonly,shared)
5235 call getenv_loc('TORPAR',torname)
5236 open (itorp,file=torname,status='old',readonly,shared)
5237 call getenv_loc('TORDPAR',tordname)
5238 open (itordp,file=tordname,status='old',readonly,shared)
5239 call getenv_loc('FOURIER',fouriername)
5240 open (ifourier,file=fouriername,status='old',readonly,shared)
5241 call getenv_loc('ELEPAR',elename)
5242 open (ielep,file=elename,status='old',readonly,shared)
5243 call getenv_loc('SIDEPAR',sidename)
5244 open (isidep,file=sidename,status='old',readonly,shared)
5246 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5247 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5248 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5249 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5250 call getenv_loc('TORPAR_NUCL',torname_nucl)
5251 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5252 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5253 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5254 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5255 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5258 #elif (defined CRAY) || (defined AIX)
5259 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5261 ! print *,"Processor",myrank," opened file 1"
5262 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5263 ! print *,"Processor",myrank," opened file 9"
5264 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5265 ! Get parameter filenames and open the parameter files.
5266 call getenv_loc('BONDPAR',bondname)
5267 open (ibond,file=bondname,status='old',action='read')
5268 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5269 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5271 ! print *,"Processor",myrank," opened file IBOND"
5272 call getenv_loc('THETPAR',thetname)
5273 open (ithep,file=thetname,status='old',action='read')
5274 ! print *,"Processor",myrank," opened file ITHEP"
5275 call getenv_loc('ROTPAR',rotname)
5276 open (irotam,file=rotname,status='old',action='read')
5277 ! print *,"Processor",myrank," opened file IROTAM"
5278 call getenv_loc('TORPAR',torname)
5279 open (itorp,file=torname,status='old',action='read')
5280 ! print *,"Processor",myrank," opened file ITORP"
5281 call getenv_loc('TORDPAR',tordname)
5282 open (itordp,file=tordname,status='old',action='read')
5283 ! print *,"Processor",myrank," opened file ITORDP"
5284 call getenv_loc('SCCORPAR',sccorname)
5285 open (isccor,file=sccorname,status='old',action='read')
5286 ! print *,"Processor",myrank," opened file ISCCOR"
5287 call getenv_loc('FOURIER',fouriername)
5288 open (ifourier,file=fouriername,status='old',action='read')
5289 ! print *,"Processor",myrank," opened file IFOURIER"
5290 call getenv_loc('ELEPAR',elename)
5291 open (ielep,file=elename,status='old',action='read')
5292 ! print *,"Processor",myrank," opened file IELEP"
5293 call getenv_loc('SIDEPAR',sidename)
5294 open (isidep,file=sidename,status='old',action='read')
5296 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5297 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5298 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5299 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5300 call getenv_loc('TORPAR_NUCL',torname_nucl)
5301 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5302 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5303 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5304 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5305 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5307 call getenv_loc('LIPTRANPAR',liptranname)
5308 open (iliptranpar,file=liptranname,status='old',action='read')
5309 call getenv_loc('TUBEPAR',tubename)
5310 open (itube,file=tubename,status='old',action='read')
5311 call getenv_loc('IONPAR',ionname)
5312 open (iion,file=ionname,status='old',action='read')
5314 ! print *,"Processor",myrank," opened file ISIDEP"
5315 ! print *,"Processor",myrank," opened parameter files"
5317 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5318 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5319 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5320 ! Get parameter filenames and open the parameter files.
5321 call getenv_loc('BONDPAR',bondname)
5322 open (ibond,file=bondname,status='old')
5323 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5324 open (ibond_nucl,file=bondname_nucl,status='old')
5326 call getenv_loc('THETPAR',thetname)
5327 open (ithep,file=thetname,status='old')
5328 call getenv_loc('ROTPAR',rotname)
5329 open (irotam,file=rotname,status='old')
5330 call getenv_loc('TORPAR',torname)
5331 open (itorp,file=torname,status='old')
5332 call getenv_loc('TORDPAR',tordname)
5333 open (itordp,file=tordname,status='old')
5334 call getenv_loc('SCCORPAR',sccorname)
5335 open (isccor,file=sccorname,status='old')
5336 call getenv_loc('FOURIER',fouriername)
5337 open (ifourier,file=fouriername,status='old')
5338 call getenv_loc('ELEPAR',elename)
5339 open (ielep,file=elename,status='old')
5340 call getenv_loc('SIDEPAR',sidename)
5341 open (isidep,file=sidename,status='old')
5343 open (ithep_nucl,file=thetname_nucl,status='old')
5344 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5345 open (irotam_nucl,file=rotname_nucl,status='old')
5346 call getenv_loc('TORPAR_NUCL',torname_nucl)
5347 open (itorp_nucl,file=torname_nucl,status='old')
5348 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5349 open (itordp_nucl,file=tordname_nucl,status='old')
5350 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5351 open (isidep_nucl,file=sidename_nucl,status='old')
5353 call getenv_loc('LIPTRANPAR',liptranname)
5354 open (iliptranpar,file=liptranname,status='old')
5355 call getenv_loc('TUBEPAR',tubename)
5356 open (itube,file=tubename,status='old')
5357 call getenv_loc('IONPAR',ionname)
5358 open (iion,file=ionname,status='old')
5359 call getenv_loc('IONPAR_NUCL',ionnuclname)
5360 open (iionnucl,file=ionnuclname,status='old')
5362 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5364 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5365 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5366 ! Get parameter filenames and open the parameter files.
5367 call getenv_loc('BONDPAR',bondname)
5368 open (ibond,file=bondname,status='old',action='read')
5369 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5370 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5371 call getenv_loc('THETPAR',thetname)
5372 open (ithep,file=thetname,status='old',action='read')
5373 call getenv_loc('ROTPAR',rotname)
5374 open (irotam,file=rotname,status='old',action='read')
5375 call getenv_loc('TORPAR',torname)
5376 open (itorp,file=torname,status='old',action='read')
5377 call getenv_loc('TORDPAR',tordname)
5378 open (itordp,file=tordname,status='old',action='read')
5379 call getenv_loc('SCCORPAR',sccorname)
5380 open (isccor,file=sccorname,status='old',action='read')
5382 call getenv_loc('THETPARPDB',thetname_pdb)
5383 print *,"thetname_pdb ",thetname_pdb
5384 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5385 print *,ithep_pdb," opened"
5387 call getenv_loc('FOURIER',fouriername)
5388 open (ifourier,file=fouriername,status='old',readonly)
5389 call getenv_loc('ELEPAR',elename)
5390 open (ielep,file=elename,status='old',readonly)
5391 call getenv_loc('SIDEPAR',sidename)
5392 open (isidep,file=sidename,status='old',readonly)
5394 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5395 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5396 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5397 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5398 call getenv_loc('TORPAR_NUCL',torname_nucl)
5399 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5400 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5401 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5402 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5403 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5404 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5405 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5406 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5407 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5408 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5409 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5410 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5411 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5414 call getenv_loc('LIPTRANPAR',liptranname)
5415 open (iliptranpar,file=liptranname,status='old',action='read')
5416 call getenv_loc('TUBEPAR',tubename)
5417 open (itube,file=tubename,status='old',action='read')
5418 call getenv_loc('IONPAR',ionname)
5419 open (iion,file=ionname,status='old',action='read')
5420 call getenv_loc('IONPAR_NUCL',ionnuclname)
5421 open (iionnucl,file=ionnuclname,status='old',action='read')
5424 call getenv_loc('ROTPARPDB',rotname_pdb)
5425 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5428 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5429 #if defined(WINIFL) || defined(WINPGI)
5430 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5431 #elif (defined CRAY) || (defined AIX)
5432 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5434 open (iscpp_nucl,file=scpname_nucl,status='old')
5436 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5441 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5442 ! Use -DOLDSCP to use hard-coded constants instead.
5444 call getenv_loc('SCPPAR',scpname)
5445 #if defined(WINIFL) || defined(WINPGI)
5446 open (iscpp,file=scpname,status='old',readonly,shared)
5447 #elif (defined CRAY) || (defined AIX)
5448 open (iscpp,file=scpname,status='old',action='read')
5450 open (iscpp,file=scpname,status='old')
5452 open (iscpp,file=scpname,status='old',action='read')
5455 call getenv_loc('PATTERN',patname)
5456 #if defined(WINIFL) || defined(WINPGI)
5457 open (icbase,file=patname,status='old',readonly,shared)
5458 #elif (defined CRAY) || (defined AIX)
5459 open (icbase,file=patname,status='old',action='read')
5461 open (icbase,file=patname,status='old')
5463 open (icbase,file=patname,status='old',action='read')
5466 ! Open output file only for CG processes
5467 ! print *,"Processor",myrank," fg_rank",fg_rank
5468 if (fg_rank.eq.0) then
5470 if (nodes.eq.1) then
5473 npos = dlog10(dfloat(nodes-1))+1
5475 if (npos.lt.3) npos=3
5476 write (liczba,'(i1)') npos
5477 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5479 write (liczba,form) me
5480 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5481 liczba(:ilen(liczba))
5482 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5484 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5486 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5487 liczba(:ilen(liczba))//'.mol2'
5488 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5489 liczba(:ilen(liczba))//'.stat'
5491 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5492 //liczba(:ilen(liczba))//'.stat')
5493 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5496 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5497 liczba(:ilen(liczba))//'.const'
5502 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5503 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5504 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5505 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5506 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5508 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5510 rest2name=prefix(:ilen(prefix))//'.rst'
5512 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5515 #if defined(AIX) || defined(PGI)
5516 if (me.eq.king .or. .not. out1file) &
5517 open(iout,file=outname,status='unknown')
5519 if (fg_rank.gt.0) then
5520 write (liczba,'(i3.3)') myrank/nfgtasks
5521 write (ll,'(bz,i3.3)') fg_rank
5522 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5527 open(igeom,file=intname,status='unknown',position='append')
5528 open(ipdb,file=pdbname,status='unknown')
5529 open(imol2,file=mol2name,status='unknown')
5530 open(istat,file=statname,status='unknown',position='append')
5532 !1out open(iout,file=outname,status='unknown')
5535 if (me.eq.king .or. .not.out1file) &
5536 open(iout,file=outname,status='unknown')
5538 if (fg_rank.gt.0) then
5539 write (liczba,'(i3.3)') myrank/nfgtasks
5540 write (ll,'(bz,i3.3)') fg_rank
5541 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5546 open(igeom,file=intname,status='unknown',access='append')
5547 open(ipdb,file=pdbname,status='unknown')
5548 open(imol2,file=mol2name,status='unknown')
5549 open(istat,file=statname,status='unknown',access='append')
5551 !1out open(iout,file=outname,status='unknown')
5554 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5555 csa_seed=prefix(:lenpre)//'.CSA.seed'
5556 csa_history=prefix(:lenpre)//'.CSA.history'
5557 csa_bank=prefix(:lenpre)//'.CSA.bank'
5558 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5559 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5560 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5561 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5562 csa_int=prefix(:lenpre)//'.int'
5563 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5564 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5565 csa_in=prefix(:lenpre)//'.CSA.in'
5566 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5569 write (iout,'(80(1h-))')
5570 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5571 write (iout,'(80(1h-))')
5572 write (iout,*) "Input file : ",&
5573 pref_orig(:ilen(pref_orig))//'.inp'
5574 write (iout,*) "Output file : ",&
5575 outname(:ilen(outname))
5577 write (iout,*) "Sidechain potential file : ",&
5578 sidename(:ilen(sidename))
5580 write (iout,*) "SCp potential file : ",&
5581 scpname(:ilen(scpname))
5583 write (iout,*) "Electrostatic potential file : ",&
5584 elename(:ilen(elename))
5585 write (iout,*) "Cumulant coefficient file : ",&
5586 fouriername(:ilen(fouriername))
5587 write (iout,*) "Torsional parameter file : ",&
5588 torname(:ilen(torname))
5589 write (iout,*) "Double torsional parameter file : ",&
5590 tordname(:ilen(tordname))
5591 write (iout,*) "SCCOR parameter file : ",&
5592 sccorname(:ilen(sccorname))
5593 write (iout,*) "Bond & inertia constant file : ",&
5594 bondname(:ilen(bondname))
5595 write (iout,*) "Bending parameter file : ",&
5596 thetname(:ilen(thetname))
5597 write (iout,*) "Rotamer parameter file : ",&
5598 rotname(:ilen(rotname))
5601 write (iout,*) "Thetpdb parameter file : ",&
5602 thetname_pdb(:ilen(thetname_pdb))
5605 write (iout,*) "Threading database : ",&
5606 patname(:ilen(patname))
5608 write (iout,*)" DIRTMP : ",&
5610 write (iout,'(80(1h-))')
5613 end subroutine openunits
5614 !-----------------------------------------------------------------------------
5617 use geometry_data, only: nres,dc
5619 ! implicit real*8 (a-h,o-z)
5620 ! include 'DIMENSIONS'
5621 ! include 'COMMON.CHAIN'
5622 ! include 'COMMON.IOUNITS'
5623 ! include 'COMMON.MD'
5626 ! real(kind=8) :: var,ene
5628 open(irest2,file=rest2name,status='unknown')
5629 read(irest2,*) totT,EK,potE,totE,t_bath
5632 ! AL 4/17/17: Now reading d_t(0,:) too
5634 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5637 ! AL 4/17/17: Now reading d_c(0,:) too
5639 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5642 read (irest2,*) iset
5646 end subroutine readrst
5647 !-----------------------------------------------------------------------------
5648 subroutine read_fragments
5652 use control_data, only:out1file
5655 ! implicit real*8 (a-h,o-z)
5656 ! include 'DIMENSIONS'
5660 ! include 'COMMON.SETUP'
5661 ! include 'COMMON.CHAIN'
5662 ! include 'COMMON.IOUNITS'
5663 ! include 'COMMON.MD'
5664 ! include 'COMMON.CONTROL'
5667 ! real(kind=8) :: var,ene
5669 read(inp,*) nset,nfrag,npair,nfrag_back
5671 !el from module energy
5672 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5673 if(.not.allocated(wfrag_back)) then
5674 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5675 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5677 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5678 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5680 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5681 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5684 if(me.eq.king.or..not.out1file) &
5685 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5686 " nfrag_back",nfrag_back
5688 read(inp,*) mset(iset)
5690 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5692 if(me.eq.king.or..not.out1file) &
5693 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5694 ifrag(2,i,iset), qinfrag(i,iset)
5697 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5699 if(me.eq.king.or..not.out1file) &
5700 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5701 ipair(2,i,iset), qinpair(i,iset)
5704 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5705 wfrag_back(3,i,iset),&
5706 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5707 if(me.eq.king.or..not.out1file) &
5708 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5709 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5713 end subroutine read_fragments
5714 !-----------------------------------------------------------------------------
5716 !-----------------------------------------------------------------------------
5720 ! implicit real*8 (a-h,o-z)
5721 ! include 'DIMENSIONS'
5722 ! include 'COMMON.CSA'
5723 ! include 'COMMON.BANK'
5724 ! include 'COMMON.IOUNITS'
5726 ! integer :: ntf,ik,iw_pdb
5727 ! real(kind=8) :: var,ene
5729 open(icsa_in,file=csa_in,status="old",err=100)
5730 read(icsa_in,*) nconf
5731 read(icsa_in,*) jstart,jend
5732 read(icsa_in,*) nstmax
5733 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5734 read(icsa_in,*) nran0,nran1,irr
5735 read(icsa_in,*) nseed
5736 read(icsa_in,*) ntotal,cut1,cut2
5737 read(icsa_in,*) estop
5738 read(icsa_in,*) icmax,irestart
5739 read(icsa_in,*) ntbankm,dele,difcut
5740 read(icsa_in,*) iref,rmscut,pnccut
5741 read(icsa_in,*) ndiff
5748 end subroutine csa_read
5749 !-----------------------------------------------------------------------------
5750 subroutine initial_write
5753 ! implicit real*8 (a-h,o-z)
5754 ! include 'DIMENSIONS'
5755 ! include 'COMMON.CSA'
5756 ! include 'COMMON.BANK'
5757 ! include 'COMMON.IOUNITS'
5759 ! integer :: ntf,ik,iw_pdb
5760 ! real(kind=8) :: var,ene
5762 open(icsa_seed,file=csa_seed,status="unknown")
5763 write(icsa_seed,*) "seed"
5765 #if defined(AIX) || defined(PGI)
5766 open(icsa_history,file=csa_history,status="unknown",&
5769 open(icsa_history,file=csa_history,status="unknown",&
5772 write(icsa_history,*) nconf
5773 write(icsa_history,*) jstart,jend
5774 write(icsa_history,*) nstmax
5775 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5776 write(icsa_history,*) nran0,nran1,irr
5777 write(icsa_history,*) nseed
5778 write(icsa_history,*) ntotal,cut1,cut2
5779 write(icsa_history,*) estop
5780 write(icsa_history,*) icmax,irestart
5781 write(icsa_history,*) ntbankm,dele,difcut
5782 write(icsa_history,*) iref,rmscut,pnccut
5783 write(icsa_history,*) ndiff
5785 write(icsa_history,*)
5788 open(icsa_bank1,file=csa_bank1,status="unknown")
5789 write(icsa_bank1,*) 0
5793 end subroutine initial_write
5794 !-----------------------------------------------------------------------------
5795 subroutine restart_write
5798 ! implicit real*8 (a-h,o-z)
5799 ! include 'DIMENSIONS'
5800 ! include 'COMMON.IOUNITS'
5801 ! include 'COMMON.CSA'
5802 ! include 'COMMON.BANK'
5804 ! integer :: ntf,ik,iw_pdb
5805 ! real(kind=8) :: var,ene
5807 #if defined(AIX) || defined(PGI)
5808 open(icsa_history,file=csa_history,position="append")
5810 open(icsa_history,file=csa_history,access="append")
5812 write(icsa_history,*)
5813 write(icsa_history,*) "This is restart"
5814 write(icsa_history,*)
5815 write(icsa_history,*) nconf
5816 write(icsa_history,*) jstart,jend
5817 write(icsa_history,*) nstmax
5818 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5819 write(icsa_history,*) nran0,nran1,irr
5820 write(icsa_history,*) nseed
5821 write(icsa_history,*) ntotal,cut1,cut2
5822 write(icsa_history,*) estop
5823 write(icsa_history,*) icmax,irestart
5824 write(icsa_history,*) ntbankm,dele,difcut
5825 write(icsa_history,*) iref,rmscut,pnccut
5826 write(icsa_history,*) ndiff
5827 write(icsa_history,*)
5828 write(icsa_history,*) "irestart is: ", irestart
5830 write(icsa_history,*)
5834 end subroutine restart_write
5835 !-----------------------------------------------------------------------------
5837 !-----------------------------------------------------------------------------
5838 subroutine write_pdb(npdb,titelloc,ee)
5840 ! implicit real*8 (a-h,o-z)
5841 ! include 'DIMENSIONS'
5842 ! include 'COMMON.IOUNITS'
5843 character(len=50) :: titelloc1
5844 character*(*) :: titelloc
5845 character(len=3) :: zahl
5846 character(len=5) :: liczba5
5848 integer :: npdb !,ilen
5852 ! real(kind=8) :: var,ene
5856 if (npdb.lt.1000) then
5857 call numstr(npdb,zahl)
5858 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5860 if (npdb.lt.10000) then
5861 write(liczba5,'(i1,i4)') 0,npdb
5863 write(liczba5,'(i5)') npdb
5865 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5867 call pdbout(ee,titelloc1,ipdb)
5870 end subroutine write_pdb
5871 !-----------------------------------------------------------------------------
5873 !-----------------------------------------------------------------------------
5874 subroutine write_thread_summary
5875 ! Thread the sequence through a database of known structures
5876 use control_data, only: refstr
5878 use energy_data, only: n_ene_comp
5880 ! implicit real*8 (a-h,o-z)
5881 ! include 'DIMENSIONS'
5883 use MPI_data !include 'COMMON.INFO'
5886 ! include 'COMMON.CONTROL'
5887 ! include 'COMMON.CHAIN'
5888 ! include 'COMMON.DBASE'
5889 ! include 'COMMON.INTERACT'
5890 ! include 'COMMON.VAR'
5891 ! include 'COMMON.THREAD'
5892 ! include 'COMMON.FFIELD'
5893 ! include 'COMMON.SBRIDGE'
5894 ! include 'COMMON.HEADER'
5895 ! include 'COMMON.NAMES'
5896 ! include 'COMMON.IOUNITS'
5897 ! include 'COMMON.TIME1'
5899 integer,dimension(maxthread) :: ip
5900 real(kind=8),dimension(0:n_ene) :: energia
5902 integer :: i,j,ii,jj,ipj,ik,kk,ist
5903 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5905 write (iout,'(30x,a/)') &
5906 ' *********** Summary threading statistics ************'
5907 write (iout,'(a)') 'Initial energies:'
5908 write (iout,'(a4,2x,a12,14a14,3a8)') &
5909 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5910 'RMSnat','NatCONT','NNCONT','RMS'
5911 ! Energy sort patterns
5916 enet=ener(n_ene-1,ip(i))
5919 if (ener(n_ene-1,ip(j)).lt.enet) then
5921 enet=ener(n_ene-1,ip(j))
5933 ist=nres_base(2,ii)+ipatt(2,i)
5935 energia(i)=ener0(kk,i)
5937 etot=ener0(n_ene_comp+1,i)
5938 rmsnat=ener0(n_ene_comp+2,i)
5939 rms=ener0(n_ene_comp+3,i)
5940 frac=ener0(n_ene_comp+4,i)
5941 frac_nn=ener0(n_ene_comp+5,i)
5944 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5945 i,str_nam(ii),ist+1,&
5946 (energia(print_order(kk)),kk=1,nprint_ene),&
5947 etot,rmsnat,frac,frac_nn,rms
5949 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5950 i,str_nam(ii),ist+1,&
5951 (energia(print_order(kk)),kk=1,nprint_ene),etot
5954 write (iout,'(//a)') 'Final energies:'
5955 write (iout,'(a4,2x,a12,17a14,3a8)') &
5956 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5957 'RMSnat','NatCONT','NNCONT','RMS'
5961 ist=nres_base(2,ii)+ipatt(2,i)
5963 energia(kk)=ener(kk,ik)
5965 etot=ener(n_ene_comp+1,i)
5966 rmsnat=ener(n_ene_comp+2,i)
5967 rms=ener(n_ene_comp+3,i)
5968 frac=ener(n_ene_comp+4,i)
5969 frac_nn=ener(n_ene_comp+5,i)
5970 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5971 i,str_nam(ii),ist+1,&
5972 (energia(print_order(kk)),kk=1,nprint_ene),&
5973 etot,rmsnat,frac,frac_nn,rms
5975 write (iout,'(/a/)') 'IEXAM array:'
5976 write (iout,'(i5)') nexcl
5978 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5980 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5981 'Max. time for threading step ',max_time_for_thread,&
5982 'Average time for threading step: ',ave_time_for_thread
5984 end subroutine write_thread_summary
5985 !-----------------------------------------------------------------------------
5986 subroutine write_stat_thread(ithread,ipattern,ist)
5988 use energy_data, only: n_ene_comp
5990 ! implicit real*8 (a-h,o-z)
5991 ! include "DIMENSIONS"
5992 ! include "COMMON.CONTROL"
5993 ! include "COMMON.IOUNITS"
5994 ! include "COMMON.THREAD"
5995 ! include "COMMON.FFIELD"
5996 ! include "COMMON.DBASE"
5997 ! include "COMMON.NAMES"
5998 real(kind=8),dimension(0:n_ene) :: energia
6000 integer :: ithread,ipattern,ist,i
6001 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6003 #if defined(AIX) || defined(PGI)
6004 open(istat,file=statname,position='append')
6006 open(istat,file=statname,access='append')
6009 energia(i)=ener(i,ithread)
6011 etot=ener(n_ene_comp+1,ithread)
6012 rmsnat=ener(n_ene_comp+2,ithread)
6013 rms=ener(n_ene_comp+3,ithread)
6014 frac=ener(n_ene_comp+4,ithread)
6015 frac_nn=ener(n_ene_comp+5,ithread)
6016 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6017 ithread,str_nam(ipattern),ist+1,&
6018 (energia(print_order(i)),i=1,nprint_ene),&
6019 etot,rmsnat,frac,frac_nn,rms
6022 end subroutine write_stat_thread
6023 !-----------------------------------------------------------------------------
6025 !-----------------------------------------------------------------------------
6026 end module io_config